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Abstract 



Understanding the collective behavior of a quantum many-body system, a system com- 
posed of a large number of interacting microscopic degrees of freedom, is a key aspect in 
many areas of contemporary physics. However, as a direct consequence of the difficultly of 
the so-called many-body problem, many exotic quantum phenomena involving extended 
systems, such as high temperature superconductivity, remain not well understood on a 
theoretical level. 

Entanglement renormalization is a recently proposed numerical method for the simulation 
of many-body systems which draws together ideas from the renormalization group and 
from the field of quantum information. By taking due care of the quantum entanglement 
of a system, entanglement renormalization has the potential to go beyond the limitations 
of previous numerical methods and to provide new insight to quantum collective phe- 
nomena. This thesis comprises a significant portion of the research development of ER 
following its initial proposal. This includes exploratory studies with ER in simple systems 
of free particles, the development of the optimisation algorithms associated to ER, and 
the early applications of ER in the study of quantum critical phenomena and frustrated 
spin systems. 

Keywords: Entanglement renormalization, quantum many-body systems, simulation 
algorithms, tensor networks, quantum entanglement, quantum information. 
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Chapter 1 



Overview 



"/t would indeed be remarkable if Nature fortified herself against further ad- 
vances in knowledge behind the analytical difficulties of the many-body prob- 
lem." 

- Max Born (1960) 

Understanding the collective behavior of quantum many-body systems, that is quantum 
systems composed of a large number of interacting microscopic degrees of freedom, is 
a central aspect in many areas of contemporary physics including particle physics and 
condensed matter physics. Solving a many-body problem to calculate, for instance, equi- 
librium properties or the outcome of a dynamic process, even while allowing for some 
degree of approximation, is often very difficult. Many techniques, both analytical and 
numeric, have been developed to address specific classes of quantum many-body prob- 
lems. Techniques based upon perturbation theory have proved incredibly successful for 
weakly interacting systems like those encountered e.g. in quantum electro-dynamics 
(QED) (MS93[ ISreOTj . For instance, perturbative calculations of electromagnetic fine 
structure constant a have been verified experimentally to astonishing accuracy; to within 
ten parts in a billion (10^^). 

In other branches of physics, notably condensed matter physics, relevant many-body sys- 
tems are often no longer in the weakly interacting regime and cannot be analysed by 
perturbation theory. Instead, without the general purpose tool of perturbation theory 
(and given that few systems of interest admit analytic solutions), much of condensed 
matter physics has been traditionally based upon a set on approximations specific to the 
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problem at hand. Unfortunately, satisfactory tools or approximations are not known for 
many models of potential interest in condensed matter theory, and these models have 
remained defiant to theoretical investigation. For instance the Hubbard model |Hub63] . 
a simple model of interacting fermions on a lattice proposed by John Hubbard in 1963, 
is thought to possibly describe high-temperature superconductivity. Despite intensive re- 
search efforts over many years, this possibility has neither been confirmed nor refuted, nor 
is even the phase diagram of the Hubbard model well understood. In general, the difficulty 
of solving many-body problems, such as the Hubbard model, has retarded progress in the- 
oretical condensed matter physics. As a consequence many exotic quantum phenomena 
involving extended systems, such as the aforementioned high-temperature superconductor 
or the spin-liquid phase of matter, remain not well understood on a theoretical level. 

Given the staggering progress made in computer technology in recent decades, numerical 
approaches to many-body problems are becoming increasingly dominant. Unfortunately, 
even with vast amounts of computational power available, large quantum many-body 
systems still cannot be addressed exactly. For instance, one would typically like to analyse 
systems with a large number of particles (often many hundreds or thousands) yet, due the 
exponential growth of the Hilbert space with particle number, only a few tens of particles 
may be analyzed exactly. As of 2009, a powerful supercomputer can analyze systems of at 
most 48 spin half particles using exact diagonalisation (ED) techniques |LL09tlRS10j . 
Despite this limitation, exact diagionalization (often combined with finite-size scaling 
techniques) remains a ubiquitous numerical tool in many areas of computational physics. 

Numerical techniques based upon Monte Carlo sampling of the multi-dimensional integrals 
which arise in quantum many-body problems have proved incredibly useful and versatile 
for studying many-body problems. These comprise a large class of algorithms collectively 
referred to as quantum Monte Carlo (QMC) |McM65tlAnd75t Cep95| . Indeed, much of the 



contemporary understanding of lattice models in D = 2 or higher dimensions has come 
from studies with QMC. Quantum Monte Carlo allows calculation of many-body effects 
in the wavefunction at the cost of statistical uncertainty that can be reduced by increasing 
the number of samples taken. While QMC techniques have proved an invaluable tool for 
many-body systems they are unfortunately only viable for a certain class of Hamiltonians. 
Specifically, QMC techniques can only be efficiently applied to study systems for which 
the Feynman path integral can be evaluated as a sum over configurations with positive 
weights. However, there is a large class of systems for which the path integral does not 
have any known representation with only positive weights and the sampling time becomes 
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exponentially large with system size; this is the famous sign problem. Two important 
classes of problems in which the sign problem is encountered are frustrated models of 
magnetism and systems of interacting fermions. Despite intensive research efforts there 
remains no general solution to the sign problem (though there have been partial solutions 
found for particular systems, for instance |HSOO[ |BHR03[ [CD04|, lAKlOj ) and many models 



of interest remain out of reach to QMC techniques. 

Another leading numerical method for many body systems is the density matrix renormal- 
ization group (DMRG) |Whi92t IWhi93] , a method based upon ideas from the real-space 
renormalization group (RG) |Fis98 ]. In the context of lattice models, the RG aims to 
obtain the physics of low energy states by grouping degrees of freedom and retaining only 
the relevant ones. Wilson's numerical renormalization group (NRG) |Wil75] provided an 
explicit prescription to implement rescaling transformations of the lattice and was suc- 
cessful in solving the Kondo problem. But it was not until the advent of White's DMRG 
algorithm in 1991 that RG methods became (and remain today) the dominant numer- 
ical approach for ID lattice systems. Typically, the DMRG algorithm can be used to 
address the ground state of large ID lattice systems (hundreds of sites) with many digits 
of precision. Unlike QMC techniques, DMRG does not suffer from the sign problem and 
can be applied to study interacting fermions and frustrated spin systems. The advent of 
the DMRG algorithm revolutionised the study of quantum systems; by allowing highly 
accurate simulation of a large class of ID lattice problems DMRG has led to a deeper 
understanding of many types of quantum systems including fermionic systems, such as 
the Hubbard model, bosonic systems, problems with impurities, and quantum dots joined 
with quantum wires. However, in D = 2 or higher dimensional quantum lattice systems, 
DMRG is limited to small system sizes (A^ 10 x 10 sites), although this is still signifi- 
cantly larger than what is possible with exact diagonalization alone. This limitation arises 
due to the failure of the matrix product state (MPS) ansatz, the tensor-network ansatz 
upon which DMRG is based, to properly capture the geometry of quantum correlations 
in higher dimensional systems. 

Though QMC and DMRG are both stunningly successful techniques in their own right, 
certain classes of many body problem remain intractable to either method; mainly frus- 
trated spin systems and interacting fermions on large 2D lattices. These systems cannot 
be analyzed with QMC due to the sign problem, nor with DMRG due to the inability 
of DMRG to analyze large 2D systems. Accordingly, there have been significant re- 
search efforts to devise new numerical many-body techniques that both (z) do not possess 
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a sign problem and (u) scale efficiently with system size in D = 2 dimensional quan- 
tum lattice systems. Several of these proposals are based upon generalizations of the 
MPS tensor network ansatz, or proposals of entirely new tensor network ansatz, that 
properly capture the correlation structure of 2D systems. Among this new breed of 
many-body techniques, based upon tensor network ansatz, are projected entangled pair 
states (PEPS) |PVC06t lJOV"'"08j . tensor entanglement renormahzation group (TERG) 
|GLW08t IGLSW09] and, the subject of this thesis, entanglement renormahzation (ER) 
|Vid07|IVid08] . 

We now proceed to outline the various aspects of entanglement renormahzation that are 
analysed in this thesis. 



1.1 Foundations of entanglement renormalization and 
the MERA 

In Chapter [2| we present a short introduction to entanglement renormalization and its 
related tensor network, the multiscale entanglement renormalization ansatz (MERA). 
The MERA is presented both from the perspective as a peculiar class of quantum circuit 
and also from the complementary perspective of the renormalization group (RG). Several 
realizations of the MERA in ID and 2D lattices are presented. The notion of the causal 
cone of the MERA is introduced and the impact on computational efficiency is discussed. 
We also discuss how spatial symmetries can be incorporated into the MERA. 



1.2 Entanglement renormalization in free fermionic 
systems 

Chapter [3] explores the performance of entanglement renormalization in systems of free 
spinless fermions on ID and 2D lattices. The implementation of ER considered makes use 
of properties of free fermions to achieve substantial reductions in computational cost, and 
also to simplify the subsequent analysis of results. The ability of ER to accurately coarse- 
grain fermionic systems in insulating, conducting and superconducting phases (whose 
corresponding ground states span all known forms of entropy scaling) is investigated. We 
examine the accuracy of the ER simulations as compared to exact results in terms of 



1.3 Entanglement renormalization in free bosonic systems 
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the truncation error of the correlation matrix, the accuracy of the ground energy and the 
accuracy of two-point correlators. The results demonstrate the ability of a coarse-graining 
transformation based upon ER to accurately capture the low-energy physics of a variety 
of free fermion systems in Z) = 1, 2 dimensional lattices. 

1.3 Entanglement renormalization in free bosonic sys- 
tems 

In Chapter |4], we investigate the performance of ER in systems of free bosons, while high- 
lighting the connection between momentum-space RG transformations and real-space RG 
transformations. The process of applying RG transforms to the system in momentum- 
space is explained and contrasted against the equivalent steps of applying real- space RG. 
A critical system and examples of relevant and irrelevant perturbations thereof are con- 
sidered. A comparison between results obtained from exact momentum-space RG and the 
results from numerical real-space RG is presented. The results of this Chapter demon- 
strate that a real-space RG transform based upon entanglement renormalization is able 
to reproduce the exact results from momentum-space RG to a high accuracy in D = 1, 2 
dimensional lattice systems, both for the critical and non-critical cases considered. Also 
demonstrated is the ability of the MERA to provide an efficient and accurate represen- 
tation of the ground state of free boson systems, thus extending to the bosonic case the 
results of Chapter [3j 

1.4 Algorithms for entanglement renormalization 

In Chapter |5] we describe how to compute expected values of local observables and two- 
point correlators from a MERA, and also present optimisation algorithms that allow the 
MERA approximation to the ground state of an arbitrary system to be obtained. A 
highlight of the algorithms is their computational cost. For an inhomogeneous lattice 
with sites, the cost scales as 0{N), whereas for translation invariant systems it drops 
to just O(logA^). Other variations of the algorithm allow us to address infinite systems, 
and scale invariant systems (e.g. quantum critical systems), at a cost independent of N. 
We also present benchmark calculations for different ID quantum lattice models, namely 
Ising, 3-level Potts, XX and Heisenberg models. We compute ground state energies. 
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magnetizations and two-point correlators throughout the phase diagram of the Ising and 
Potts models, which includes a second order quantum phase transition. We find that, 
at the critical point of an infinite system, the error in the ground state energy decays 
exponentially with the refinement parameter whereas the two-point correlators remain 
accurate even at distances of millions of lattice sites. We then extract critical exponents 
from the order parameter and from two-point correlators. Finally, we also compute a 
MERA that includes the first excited state, from which the energy gap can be obtained 
and seen to vanish as at criticality. 

1.5 Entanglement renormalization and quantum crit- 
icality 

Chapter [6] explains how to use the MERA to investigate quantum critical systems that 
are invariant under changes of scale. An algorithm is presented which, given a critical 
Hamiltonian, details how to compute a scale invariant MERA for its ground state. Then, 
starting from the scale invariant MERA, a procedure is described to identify the scaling 
operators/dimensions of the theory. A closed expression for two-point and three-point 
correlators is derived, and a connection between the MERA and conformal field theory 
(which can be used to readily identify the continuum limit of a critical lattice model) is 
established. Finally, benchmark calculations for the Ising and Potts models are presented. 

1.6 Entanglement renormalization in two spatial di- 
mensions 

In Chapter [7] we build upon the MERA schemes introduced in Chapter [2} and also the 
optimisation algorithms presented in Chapter [5} to describe an implementation of the 
MERA that allows us to consider 2D systems of arbitrary size, including infinite systems. 
In this way we demonstrate the scalability of entanglement renormalization in two spatial 
dimensions and decisively contribute to establishing the MERA as a competitive approach 
to systematically address 2D lattice models. We also demonstrate the performance of the 
scheme by analysing the 2D quantum Ising model, for which we obtain accurate estimates 
of the ground state energy and magnetizations, as well as two-point correlators (shown to 
scale polynomially at criticality), the energy gap, and the critical magnetic field and beta 
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exponent. Finally, we discuss how the use of disentanglers affects the simulation costs, by 
comparing the MERA with a tree tensor network (TTN). 

1.7 The spin- 1 Kagome lattice Heisenberg model with 
entanglement renormalization 

Chapter [8] examines a long-standing problem in condensed matter physics, the nature of 
the ground-state of the spin-| kagome lattice Heisenberg model (KLHM), using entangle- 
ment renormalization. Progress in understanding the KLHM has been hindered, as with 
many other models of frustrated antiferromagnets, by the inapplicability of quantum 
Monte Carlo methods due to the sign problem. This investigation is the first demon- 
stration of the utility of entanglement renormalization to study 2D lattice models that 
are beyond the reach of quantum Monte Carlo techniques. After describing a scheme 
for ER on the Kagome lattice with periodic boundary conditions, we address lattices of 
= {36, 144, oo} sites. We then analyse bond energies, two-point and bond-bond corre- 
lators of the ground state on the various lattice sizes. The results give strong numerical 
evidence in favor of a VBC ground state for the KLHM. 
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Chapter 2 

Foundations of entanglement 
renormalization and the MERA 



2.1 Introduction 

Entanglement renormalization (ER) |Vid07] is a numerical technique based on locally 
reorganizing the Hilbert space of a quantum many-body system with the aim to reduce 
the amount of entanglement in its wave function. It was introduced to address a major 
computational obstacle in real space renormalization group (RG) methods |Wil75[ rWhi92l 



IWhi93t IMW941 IMW96j . responsible for limitations in their performance and range of 
applicability, namely the proliferation of degrees of freedom that occurs under successive 
applications of a RG transformation. 

Entanglement renormalization is built around the assumption that, as a result of the local 
character of physical interactions, some of the relevant degrees of freedom in the ground 
state of a many-body system can be decoupled from the rest by unitarily transforming 
small regions of space. Accordingly, unitary transformations known as disentanglers are 
applied locally to the system in order to identify and decouple such degrees of freedom, 
which are then safely removed and therefore do no longer appear in any subsequent coarse- 
grained description. This prevents the harmful accumulation of degrees of freedom and 
thus leads to a sustainable real space RG transformation, able to explore arbitrarily large 
ID and 2D lattice systems, even at a quantum critical point. It also leads to the multi- 
scale entanglement renormalization ansatz (MERA), a variational ansatz for many-body 
states |Vid08j . 
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The MERA, based in turn on a class of quantum circuits, is particularly successful at 
describing ground states at quantum criticality |Vid07l IVidOSl lEVlObl lEVlOal ICDR081 
l(;MF08t IPEVOQi IEV09bj or with topological order |AV08[ IKRVOQj . From the computa- 
tional viewpoint, the key property of the MERA is that it can be manipulated efficiently, 
due to the causal structure of the underlying quantum circuit |Vid08j . As a result, it is 
possible to efficiently evaluate the expected value of local observables, and to efficiently 
optimize its tensors. Thus, well-established simulation techniques for matrix product 
states, such as energy minimization |VPC04l IPVCOGj or simulation of time evolution 
[VidOai IVid04j . can be readily generahzed to the MERA |DE()08[ IRMV08j . 

The goal of this Chapter is to provide a short introduction to entanglement renormaliza- 
tion and the MERA, establishing the notation and terminology that shall be used in the 



rest of this thesis. Specifically, Sect. 2.3 introduces the MERA from the perspective of 



quantum circuits while Sect. 2. 4| offers a complementary interpretation of the MERA in 



terms of the renormalization group (RG). Several realizations of the MERA in ID and 
2D lattices are discussed in Sect. 2]5]and Sect. 2^ details how spatial symmetries can be 
incorporated into the MERA. 



2.2 The MERA 



Let C denote a D-dimensional lattice made of sites, where each site is described by a 
Hilbert space V of finite dimension d, so that Yc — V®^. The MERA is an ansatz to 
describe certain pure states |^) G of the lattice or, more generally, subspaces Yu C Yc- 

There are two useful ways of thinking about the MERA that can be used to motivate its 
specific structure as a tensor network, and also help understand its properties and how 
the algorithms ultimately work. One way is to regard the MERA as a quantum circuit 
C whose output wires correspond to the sites of the lattice £ |Vid08] . Alternatively, 
we can think of the MERA as defining a coarse-graining transformation that maps C 
into a sequence of increasingly coarser lattices, thus leading to a renormalization group 
transformation |Vid07] . Next we briefiy review these two complementary interpretations. 



2.2 The MERA 
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Quantum Circuit: 





Figure 2.1: Quantum circuit C corresponding to a specific realization of the MERA, 



namely the binary ID MERA of Fig. |2.2[ In this particular example, circuit C is made 
of gates involving two incoming wires and two outgoing wires, p = = Pout = 2. Some 
of the unitary gates in this circuit have one incoming wire in the fixed state |0) and can 
be replaced with an isometry w of type (1,2). By making this replacement, we obtain the 



isometric circuit of Fig. 2.2 
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Binary ID MERA: 

k=J,2,...,XT 




h '2 '3 "4 'i "« '7 I9 '10 '11 '12 '13 '14 '15 'l« 




Figure 2.2: Top: Example of a binary ID MERA for a lattice C with iV = 16 sites. It 
contains two types of isometric tensors, organized in T = 4 layers. The input (output) 
wires of a tensor are those that enter it from the top (leave it from the bottom). The 
top tensor is of type (1, 2) and the rank xt of its upper index determines the dimension 
of the subspace Yjj C represented by the MERA. The isometries w are of type (1, 2) 
and are used to replace each block of two sites with a single effective site. Finally, the 
disentanglers u are of type (2,2) and are used to disentangle the blocks of sites before 
coarse-graining. Bottom: Under the renormalization group transformation induced by 
the binary ID MERA, three-site operators are mapped into three-site operators. 



2.3 Quantum circuit 
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2.3 Quantum circuit 

As a quantum circuit C, the MERA for a pure state |\E') G Yc is made of quantum 
wires, each one described by a Hilbert space V, and unitary gates u that transform the 



unentangled state lO)*^^ into (see Fig. 2.1). 



In a generic case, each unitary gate u in the circuit C involves some small number p of 
wires, 

u:Y'^P^Y^P, u^u = uu^ = I, (2.1) 

where I is the identity operator in V®^. For some gates, however, one or several of the 
input wires are in a fixed state |0). In this case we can replace the unitary gate u with 
an isometric gate w 

W.Yin^Yout, w^w = h,„, (2.2) 

where Vj„ = V®^'" is the space of the Pin input wires that are not in a fixed state |0) and 
Yout — V^^""' is the space of the Pout = P output wires. We refer to w as a {pin,Pout) gate 
or tensor. 



Fig. 2.2 shows an example of a MERA for a ID lattice C made of A^ = 16 sites. Its tensors 



are of types (1, 2) and (2, 2). We call the (1, 2) tensors isometrics w and the (2, 2) tensors 



disentanglers u for reasons that will be explained shortly, and refer to Fig. |2.2| as a binary 
ID MERA, since it becomes a binary tree when we remove the disentanglers. Most of the 
early work for ID lattices |Vid07l IVidOSl lEVlObi lEVlOal IDEO081 IRMVOSj has been done 
using the binary ID MERA. However, there are many other possible choices. In Chapters 



^ and ml for instance, we will mostly use the ternary ID MERA of Fig. 2.3, where the 



isometrics w are of type (1, 3) and the disentanglers u remain of type (2, 2). Fig. 2.4 makes 



more explicit the meaning of Eq. |2.2| for these tensors. Notice that describing tensors 
and their manipulations by means of diagrams is fully equivalent to using equations and 
often much more clear. 



Eq. 2.2 encapsulates a distinctive property of the MERA as a tensor network: each of its 



tensors is isometric (notice that Eq. 2.1 is a particular case of Eq. 2.2). A second key 
feature of the MERA refers to its causal structure. We define the past causal cone of an 
outgoing wire of circuit C as the set of wires and gates that can affect the state on that 
wire. A quantum circuit C leads to a MERA only when the causal cone of an outgoing 
wire involves just a constant (that is, independent of A^) number of wires at any fixed past 
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Figure 2.3: Top: Example of ternary ID MERA (rank xt, T 
sites. It differs from the binary ID MERA of Fig 

(1,3), so that blocks of three sites are replaced with one effective site. 



3) for a lattice of 18 



2.2| in that the isometries are of type 

Bottom: As a 



result, two-site operators are mapped into two-site operators during the coarse-graining. 



time (Fig. 2.5). We refer to this property by saying that the causal cone has a bounded 
'width'. 

The usefulness of the quantum circuit interpretation of the MERA will become clear in 
Chap. |5| in the context of computing expected values for local observables. There, the 



two defining properties, namely Eq. 2.2 and the peculiar structure of the causal cones of 



C, will be the key to making such calculations efficient. 



2.4 Renormalization group transformation 

Let us now review how the MERA defines a coarse-graining transformation for lattice 
systems that leads to a real-space renormalization group scheme, known as entanglement 
renormalization |Vid07j . 



2.4 Renormalization group transformation 
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U 
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P 



IV 



\ lap \ h'^' m 



Figure 2.4: The tensors which comprise a MERA are constrained to be isometric, cf. 
Eq. 2.2[ The constraints for the isometrics w and disentanglers u of the ternary MERA 
can be equivalently expressed (i) diagramatically or (ii) with equations. In this thesis we 
will mostly use the diagramatic notation, which remains simple for complicated tensor 
networks. 



We start by grouping the tensors in the MERA into T log different layers, where each 
layer contains a row of isometrics w and a row of disentanglers u. We label these layers 
with an integer r = 1, 2, ■ ■ ■ T, with r = 1 for the lowest layer and with increasing values 
of r as we climb up the tensor network, and denote by Ur the isometric transformation 
implemented by all tensors in layer r, see Figs. 2.2| and |2.3 Notice that the incoming 



wires of each Ur define the Hilbert space of a lattice C-r with a number of sites that 
decreases exponentially with r (specifically, as N2~'^ and A^S""^ for the binary and ternary 
ID MERA). That is, the MERA implicitly defines a sequence of lattices 

£o ^ ^1 ^ ^ ^T, (2.3) 

where Cq = C\s the original lattice, and where we can think of lattice Cr as the result of 
coarse-graining lattice £r-i- 



Specifically, as illustrated in Figs. 2^ and 2^ this coarse-graining transformation is 
implemented by the operator Ul that maps pure states of the lattice £r-i into pure 
states of the lattice 

Ul:\c.^-,^Ncr. (2.4) 
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Figure 2.5: The past causal cone of a group of sites in £o = is the subset of wires and 
gates that can affect the state of those sites. The example shows the causal cone of a pair 
of nearest neighbor sites of Cq for the ternary ID MERA. Notice that for each lattice 
r = 0, 1, 2, 3, 4, the causal cone involves at most 2 sites. This can be seen to be the case 
for any pair of contiguous sites of Cq. We refer to this property by saying that the causal 
cones of the MERA have bounded width. 



and that proceeds in two steps (Fig. 2.6). Let us partition the lattice Cr-i into blocks of 
neighboring sites. The first step consists of applying the disentanglers u on the boundaries 
of the blocks, aiming to reduce the amount of short range entanglement in the system. 
Once (part of) the short-range entanglement between neighboring blocks has been re- 
moved, the isometrics w are used in the second step to map each block of sites of lattice 
Cr-i into a single effective site of lattice 

By composition, we obtain a sequence of increasingly coarse-grained states, 

|^o)^|^i)^ ^I^t), (2.5) 

for the lattices {Cq, Ci^ ■ ■ ■ ,£t}, where \^r) = U]. \^t-i) and |\I'o) = |^) is the original 
state. Overall the MERA corresponds to the transformation U = U1U2 ■ ■ ■ Ut, 

U:Yc,^Yco, (2.6) 

with I^^Jq) = U I^t)- 

Regarding the MERA from the perspective of the renormalization group is quite instruc- 
tive. It tells us that this ansatz is likely to describe states with a specific structure of 
internal correlations, namely, states in which the entanglement is organized in different 
length scales. Let us briefly explain what we mean by this. 



2.4 Renormalization group transformation 
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(i) Binary M ERA scheme for ID lattice 
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(ii) Ternary MERA scheme for ID lattice 
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Figure 2.6: Detailed description of the real-space renormalization group transformation 
for ID lattices induced by {%) the binary ID MERA and {ii) the ternary ID MERA. 
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(i) 2x2 MERA scheme for 2D square lattice 
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(ii) 3x3 MERA scheme for 2D square lattice 

Original Lattice Apply Disentanglers 



Coarse-grained Lattice 



oai 



. EE) • gj^ 



r Apply Isometries 



















• 






• 


•> 






o 




O 


o 




o 


• 


















If 




» 






e 


















> 

















Figure 2.7: Detailed description of the real-space renormalization group transformation 
for a 2D square lattice induced by two possible realizations of the MERA, generalizing 



the ID schemes of Fig. |2.6[ In the first case the isometries map a block of 2 x 2 sites 
into a single site, which can be seen to imply that the natural size of a local operator, 
equivalently the causal width of the scheme, is 3 x 3 sites. In the second case the isometries 
map a block of 3 x 3 sites into a single site and the natural size of a local operator is 
2x2. As a result, the computational cost in the second scheme is much smaller than in 
the first scheme. 
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We say that the state contains entanglement at a given length scale A if by applying 
a unitary operation (i.e. a disentangler) on a region R of linear size A, we are able to 
decouple (i.e. disentangle) some of the local degrees of freedom, that is, if we are able to 
convert the state into a product state ® |0), where |0) is the state of the local 
degrees of freedom that have been decoupled and is the state of the rest of the system. 
[Here we assumed, of course, that the decoupling is not possible with a unitary operation 
that acts on a subregion R' of the region R, where the size A' of R' is smaller than the 
size of R, A' < A] . 

What makes the MERA useful is that the entanglement in most ground states of local 
Hamiltonians seems to decompose into moderate contributions corresponding to different 
length scales. We can identify two behaviors, depending on whether the system is in 
a phase characterized by symmetry-breaking order or by topological order (see |LW05j 
and references therein). In systems with symmetry-breaking order, ground-state entan- 
glement spans all length scales A smaller than the correlation length ^ in the system — 
and, consequently, at a quantum critical point, where the correlation length ^ diverges, 
entanglement is present at all length scales |Vid07] . In a system with topological order, 
instead, the ground state displays some form of (topological) entanglement affecting all 
length scales even when the correlation length vanishes |AV08l IKRV09] . 



2.5 Choose your MERA 

We have introduced the MERA as a tensor network originating in a quantum circuit. Its 
tensors have incoming and outgoing wires/indices according to a well-defined direction of 
time in the circuit. Therefore, a MERA can be regarded as a tensor network equipped 
with a (fictitious) time direction and with two properties: 



Its tensors are isometric (Eq. 2.2) 



Past causal cones have bounded width (Fig. 2.5) 



From a computational perspective, these are the only properties that we need to retain. 
In particular, there is no need to keep the vector space dimension of the quantum wires 
(equivalently, of the sites in the coarse-grained lattice) constant throughout the tensor 
network. Accordingly, we will consider a MERA where the vector space dimension of a 
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site of lattice denoted Xr, may depend on the layer r (this dimension could also be 
different for each individual site of layer r, but for simplicity we will not consider this 
case here). Notice that xo = d corresponds to the sites of the original lattice C 

Bond dimension. — Often, however, the sites in most layers will have the same vector 
space dimension (except, for instance, the sites of the original lattice £, with xo = d, 
or the single site of the top lattice Ct, with dimension xt)- In this case we denote the 
dominant dimension simply by and we refer to the MERA as having bond dimension 
X- The computational cost of the algorithms described in subsequent sections is often 
expressed as a power of the bond dimension x- 

Rank. — We refer to the dimension xt of the space V^^ (corresponding to the single 
site of the uppermost lattice Ct) as the rank of the MERA. For xt = ^, the MERA 
represents a pure state |\E') G V^. More generally, a rank xt MERA encodes a xt- 
dimensional subspace Yu C V^. For instance, given a Hamiltonian H on the lattice 
£, we could use a rank xt MERA to describe the ground subspace of H (assuming it 
had dimension Xt)] or the ground state of H (if it was not degenerate) and its xt ~ 1 
excitations with lowest energy. The isometric transformation U in Eq. |2.6 can be used 
to build a projector P = UU\ 

p.,Yc^Yc, P' = P, tr(P) = xr, (2.7) 
onto the subspace Yu ^ V^. 

Given the above definition of the MERA, many different realizations are possible depend- 
ing on what kind of isometric tensors are used and how they are interconnected. We have 
already met two examples for a ID lattice, based on a binary and ternary underlying tree. 



Fig. 2/7 shows two schemes for a 2D square lattice. It is natural to ask, given a lattice 
geometry, what realization of the MERA is the most convenient from a computational 
point of view. A definitive answer to this question does not seem simple. An important 
factor, however, is given by the fixed-point size of the support of local observables under 
successive RG transformations — which corresponds to the width of the past causal cones. 

Support of local observables. — In each MERA scheme, under successive coarse- 
graining transformations a local operator eventually becomes supported in a characteristic 
number of sites. This is the result of two competing effects: disentanglers u tend to 
extend the support of the local observable (by adding new sites at its boundary), whereas 
the isometries w tend to reduce it (by transforming blocks of sites into single sites). For 
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instance, in the binary ID MERA, local observables end up supported in three contiguous 



sites (Fig. 2.2), whereas in the ternary ID MERA local observables become supported in 



two contiguous sites (Fig. 2.3). 



Therefore, an important difference between the binary and ternary ID schemes is in 
the natural support of local observables. This can be seen to imply that the cost of a 
computation scales as a larger power of the bond dimension x for the binary scheme than 
for the ternary scheme, namely as 0{x^) compared to 0{x^). However, it turns out that 
the binary scheme is more effective at removing entanglement, and as a result a smaller 
X is already sufficient in order to achieve the same degree of accuracy in the computation 
of, say, a ground state energy. In the end, we find that for the ID systems analyzed in 



Sect. |5.5[ the two effects compensate and the cost required in both schemes in order to 
achieve the same accuracy is comparable. On the other hand, in the ternary ID MERA, 
two-point correlators between selected sites can be computed at a cost 0{x^), whereas 
analogous calculations in the binary ID MERA are much more expensive. Therefore in 
any context where the calculation of two-point correlators is important, the ternary ID 
MERA is a better choice. 

The number of possible realizations of the MERA for 2D lattices is greater than for ID 



lattices. For a square lattice, the two schemes of Fig. 2.7 are obvious generalizations of 
the above ones for ID lattices. The first scheme, proposed in |EV10b] (see also |CDR08] ). 
involves isometrics of type (1,4) and the natural support of local observables is a block of 
3x3 sites. The second scheme involves isometrics of type (1,9) and local observables end 
up supported in blocks of 2 x 2 sites. Here, the much narrower causal cones of the second 
scheme leads to a much better scaling of the computational cost with x, only 0{x^^) 
compared to 0(x^^) for the first scheme. 

Another remark in relation to possible realizations concerns the type of tensors we use. 
So far we have insisted in distinguishing between disentanglers u (unitary tensors of type 
p p) and isometrics w (isometric tensors of type 1 p'). We will continue to use 
this terminology throughout this paper, but we emphasize that a more general form of 
isometric tensor, e.g. of type (2,4), that both disentangles the system and coarse-grains 
sites, is possible and is actually used in some realizations |EV09b] . 
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(i) Non tra nslationa I ly invariant MERA (ii)Translationally invariant IV1 ERA 




(iii) Scale invariant IVIERA 




Figure 2.8: Ternary ID MERA in the presence of space symmetries, (i) In order to repre- 
sent an inhomogeneous state/subspace, all disentanglers u and isometries w are different 
(denoted by different colouring). Notice that there are disentanglers (isometries) in 
the first layer, N/9 in the second, and more generally N/3'^ in layer r, so that the total 
number of tensors is 2A^ 1/3^ < 2iV. Therefore the total number of parameters 

required to specify the MERA is proportional to the size N of the lattice C (ii) In order 
to represent a state/subspace that is invariant under translations, we choose all disen- 
tanglers and isometries on a given layer of the MERA to be the same. In this case the 
MERA is completely specified by O(logA^) disentanglers and isometries. (iii) In a scale 
invariant MERA, the same disentangler and isometry is in addition used in all layers. 



2.6 Exploiting symmetries 

Symmetries have a direct impact on the efficiency of computations, because they can be 
used to drastically reduce the number of parameters in the MERA. Important examples 



are given by space symmetries, such as translation and scale invariance, see Fig. 2.8 



The MERA is made of 0{N) disentanglers and isometries. In order to describe an in- 
homogeneous state |\I') G Yc or subspace Yu ^ V^, all these tensors are chosen to be 
different. Therefore, for fixed x fhe number of parameters in the MERA scales linearly 
in N. 



However, in the presence of translation invariance, one can use a translation invariant 
MERA, where we choose all the disentanglers u and isometries w of any given layer r to 
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be the same, thus reducing the number of parameters to 0(log A^) (if there are T ^ log 
layers). We emphasize that a translation invariant MERA, as just defined, does not 
necessarily represent a translation invariant state |\E') G Yc or subspace Yu ^ V^. The 
reason is that different sites of C are placed in inequivalent positions with respect to 
the MERA. As a result, often the MERA can only approximately reproduce translation 
invariant states/subspaces, although the departure from translation invariance is seen to 
typically decrease fast with increasing x- In order to further mitigate inhomogeneities, we 
often consider an average of local observables/reduced density matrices over all possible 
sites, as will be discussed in Chap. |5} 

In systems that are invariant under changes of scale, we will use a scale invariant MERA, 
where all the disentanglers and isometrics can be chosen to be the same and we only need 
to store a constant number of parameters. The scale invariant MERA is useful to represent 
the ground state of some quantum critical systems |Vid07] and the ground subspace of 
systems with topological order at the infrared limit of the RG flow |AV08t IKRVOQj . 

A reduction in parameters (as a function of x) is also possible in the presence of internal 
symmetries, such as U{1) (e.g. particle conservation) or SU{2) (e.g. spin isotropy). 
Exploitation of internal symmetries is discussed in Ref. |SPV09] and shall not be further 
considered in this thesis. 



2.7 Conclusion 

In this Chapter, the MERA has been introduced both from the perspective of a perculiar 
class of quantum circuit and from the perspective of an RG transformation of the lattice. 
Several implementations of the MERA for D = 1 and D = 2 systems have been described 
and we have briefly discussed how the presence of symmetries can be exploited to reduce 



the computational cost of the MERA. The binary MERA of Fig. 2.2 and the 4-to-l 



MERA scheme of Fig. 2/7 shall be used in Chapters [3] and |4] to explore the performance 



of the MERA in D = 1, 2 dimensional lattices of free fermions and free bosons. 
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Chapter 3 



Entanglement renormalization in free 
fermionic systems 



3.1 Introduction 

The renormalization group (RG), concerned with the change of physics with the obser- 
vation scale, is among the main ideas underlying the theoretical structure of statistical 
mechanics and quantum field theory, and is of central importance in the modern formula- 
tion of critical phenomena and phase transitions |Fis98] . Its influence extends well beyond 
the conceptual domain: RG transformations are also the basis of numerical approaches 
to the study of strongly correlated many-body systems. 

In a lattice model, a real-space RG transformation produces a coarse-grained system by 
first joining the lattice sites into blocks and then replacing each block with an effective 
site |KGH+67j . Two very natural requirements for a RG transformation are: (i) it should 
preserve the long-distance physics of the system; (m) when this physics is invariant under 
changes of scale, the system should be a fixed point of the RG transformation. 

For the important case of a quantum system at zero temperature, the first requirement is 
fulfilled if, as determined by White in his density matrix renormalization group (DMRG) 
[Whi92l IWhi93j . the vector space of the effective site retains the local support of the 
ground state. Entanglement renormalization (ER) has been proposed in order to simul- 
taneously meet the second requirement. By using disentanglers, ER aims to produce a 
coarse-grained lattice locally identical to the original one, in the sense that their sites have 
the same vector space dimension. When this is accomplished, the original system and its 
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coarse-grained version can be meaningfully compared, e.g. through their Hamiltonians or 
ground state properties, leading to a proper real-space RG flow. 

In this Chapter we explore the performance of ER in systems of free spinless fermions 
on ID and 2D lattices using the binary MERA and 4-to-l MERA schemes described in 
Figs. 2^ and 2/7 of Chap. |2} Speciflcally, we consider systems specifled by the quadratic 



Hamiltonian 

H = ^[ajfls + o^ar — 'j{alal + a^Or)] — 2A ^ alar, (3.1) 

(rs) r 

where A and 7 are the chemical and pairing potentials and the flrst sum involves only 
nearest neighbors. In spite of its simplicity, Hamiltonian H contains a rich phase diagram 
as a function of A and 7, including insulating, conducting and superconducting phases 
[LDY"'"06] . Importantly, the corresponding ground states span all known forms of entropy 
scaling |LDY"'"06l IBCS06] . In addition, H can be diagonalized through linear (Fourier and 
Bogoliubov) transformations of the fermion operators a and while, by Wick's theorem, 
all properties of its Gaussian ground state I^E'gs) can be extracted from the two-point 
correlators (ajcis) and (arCLs). Then, provided that our RG transformation also maps 
fermion modes linearly, the entire analysis can be conducted in the space of two-point 
correlators and quadratic Hamiltonians of fermionic modes, as represented by A^ x A^ 



matrices. Hence quadratic fermionic models such as Eq. 3.1) offer an appealing testing 
ground for ER, one where computational costs have been greatly simplifled (e.g. H can 
be diagonalized exactly with just 0{N^) operations) while keeping a rich variety of non- 
trivial ground state structures. 



3.2 ER applied to free fermions 

We start by rephrasing, in the language of correlation matrices, the process of coarse- 
graining a D-dimensional (hypercubic) lattice. We assume that the system is in the 
ground state |\E'gs) of H, which we compute using standard analytic techniques (see e.g. 
[LDY"'"06] ). It is convenient to redraw the hypercubic lattice so that each site contains 
P = fermion modes for some integer p. Then a hypercube of 2^ sites deflnes a block 
that contains P2^ modes. The goal of the RG transformation is to replace this block 
with just one effective site made of P' modes, with P' < P2^. We would like to have 
P' = P, so that the sites of the coarse-grained and original lattices are identical and we 
can compare the corresponding Hamiltonians or ground-state reduced density matrices. 
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However, in the coarse-graining step only modes of the block that are disentangled from 
the rest of the system can be removed. As a result, P' often must be larger than P. 

For the sake of simplicity, we continue the analysis for the case of a ID lattice. Let us 



temporarily replace the N spinless fermion operators d in Eq. 3.1 with 2N (self-adjoint) 
Majorana fermion operators c, 

C2r-l=Cir + al, C2r = -^—^ ■ (3.2) 

% 

The ground state \^Ga) is then completely specified by 

where F, henceforth referred to as the correlation matrix, is real and antisymmetric. 
Similarly, the reduced density matrix pcs ioT a block made of 2 sites, that is with L = 2P 
spinless modes (equivalently, 2L Majorana modes) is described by a 2L x 2L submatrix 
Tl of F. This matrix is brought into (block) diagonal form by a special orthogonal 
transformation V, 



r=l 



Vr 
-Vr 



V e S0(2L), (3.4) 



where < f < 1 are the eigenvalues of F^, each one associated with a pair of Majorana 
fermions. These pairs recombine into L spinless fermions in a product state |VLRK03l 
ILRV04j 



L 



1+Vr 







PGs = 00f^. = (X)| Q 1^ h (3.5) 



r=l r=l 



where Qr, the state of a spinless fermion mode, is mixed if f ^ < 1 and pure if Vr = 1. 
Notice that since the ground state |\1/gs) is a pure state, a mode in a mixed state must be 
entangled with modes outside the block, whereas a mode in a pure state is unentangled 
from the rest of the system. We build an effective site by removing from the block, or 
projecting out from F^, all the modes that are unentangled (pure), and just keeping those 
P' modes that are entangled (mixed). In this way, the coarse-grained lattice retains the 
ground state properties. 
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The equivalence of Hilbert space truncation and the removal of 

fermionic modes 

Here we describe the process of coarse-graining a lattice by replacing blocks of sites 
with effective sites. We show that the truncation of the Hilbert space of a given block 
can be implemented by eliminating some of the modes in that block. We consider 
a fermionic lattice system in its (gaussian) ground state |\I/gs)- Let V be the vector 
space of a block containing L modes and let pcs denote the reduced density matrix 
of |\E')(3g on the block. We assume that the support of is concentrated in 

a subpace V C V. Then, following White |Whi92t IWhi93] . the optimal coarse- 
graining of the block is obtained by defining an effective site s' with vector space 
= v. In our case, pos is the tensor product of density matrices Qr for individual 
modes |VLRK03[ ILRV04] . 



Pos = 6^Qr = 6^\ I . (3.6) 




Suppose that the first P' modes are in a mixed state and the remaining L — P' 
modes are in a pure state. Then we can write 

P' L 

PgS = (0 ^r) ® ( (g) ^r) = ® TT, (3.7) 

r=l r=P'+l 

where a is a mixed state with rank 2^' whereas vr is a projector with rank 1. Let V = 
V°" ® be a tensor factorization of V such that cr = trvT(p). The key observation 
is that = V, and that p and a have the same none-vanishing eigenvalues. 
Therefore, we have two equivalent ways of constructing the space V'^' for the effective 
site s' while preserving the support of the ground state density matrix Pqq. On the 
one hand, V^' can be obtained by projecting V on the support V of p. On the 
other, V can also be build by factorizing the space V into two factor spaces V"" and 
V^, and by then tracing out the second factor, corresponding to modes in a pure 
state. Both constructions lead to an equivalent effective lattice. Finally, tracing out 
the factor space of mode r corresponds, in the language of correlation matrices Tl, 



to removing the rth row and a rth column of VTlV"^ in Eq. 3^, process to which 
we referred to as projecting out the mode. 
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Figure 3.1: Top: A block of two sites (four modes) is coarse-grained into an effective 
site by first applying disentanglers U across the boundary of the block and then using 
isometry W to project out two modes. Bottom: Same RG transformation written in the 
language of correlation matrices, Eq. |3.9 



The key idea of ER, see Fig. 3.1, is to use disentangling unitary transformations, or 
disentanglers, to diminish P' by increasing the number of modes in the block that are 
unentangled from the rest of the system. A disentangler is implemented through a special 
orthogonal matrix U G S0(2L) that acts on two neigboring sites across the boundary of 
the block, wheareas the coarse-graining is implemented by an isometry W = RYpi that 
selects the P' spinless fermion modes to be kept in the effective site, where R G S0(2L) 
and 



e 



gr 




' 1 r <P' 


, 9r= ' 


^ r>P' 


-gr _ 



(3.8) 



Let r describe three consecutive blocks. Then the correlation matrix F'r for the effective 



site reads (Fig. 3.1 ) 



r'^ = w^{u® u)^ FsL {u © u) w. 



(3.9) 
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Similarly, the correlation matrix F'^ for the modes to be removed is 

T'l = W^ {U®U)^r3LiU®U)W, (3.10) 

1 



W = R{Yl-Yp,), Yl = 



-1 



(3.11) 



Our goal is to maximize the purity of the modes to be projected out, so that they be- 
come as unentangled as possible. The sum of their purities, Yl,^=p'+i'^r-, is half of the 
antisymmetric trace of f'^, iiiV'iYl). Consequently, U and W are obtained from the 
optimization 

max ti(v'rY}\, (3.12) 
that we address through a sequence of alternating optimizations for U and R |EV09aj . 

Then, given the correlation matrix F for I^E^gs)? the RG transformation is implemented 
in three steps: (i) first a submatrix T^l for three consecutive blocks is extracted from 



F; (u) then disentangler U and isometry W are computed using the optimization (3.12) 
while keeping P' = P modes in the effective site; {iii) finally, U and W are used to 
transform the original iV-mode system into a coarse-grained system with just N/2 modes 
and effective correlation matrix T^^\ Some of the modes that are removed are still slightly 
mixed. Their mixness = 1 — f r quantifies the errors introduced. Iteration of the RG 
transformation produces a sequence of increasingly coarse-grained lattices, described by 
correlation matrices {F'-^-', F^^\ ■ ■ ■ }. The corresponding disentanglers {f/^^', ^Z*-^-*, • ■ " } 
and isometrics W^'^\ ■ ■ ■ } constitute the multi-scale entanglement renormalization 

ansatz (MERA) |Vid08] for the ground state I^^gs)- 



3.3 Results and discussion 



We have applied the present RG approach to Hamiltonian (3.1) in the thermodynamic 
limit N ^ oo. First we consider ID systems, where the whole (7, A) plane can be mapped 
into the quantum spin XY model using a Jordan- Wigner transformation, (i) For the line 
7 = 1 [equivalent to the quantum spin Ising model] we consider P = 2 modes per site 
and apply 13 iterations of the RG transformation, so that a final effective site (with just 
P = 2 modes) corresponds to 2 x 2^^ = 16384 modes of the original system. At the 
critical point A = 1, which is the most demanding, the mixedness of the removed modes 
is at most er = 1.2 x 10~^. The effect on local observables, even after the 13 iterations, 



3.3 Results and discussion 



31 



Ising Model XX Model 
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Figure 3.2: Scaling of the entanglement entropy 5*^ |VLRK03t ILRV04j in ID systems. 
Left: Quantum Ising model, 7 = 1. Bold (solid/dotted) lines represent entanglement at 
criticality, A = 1. The system is an entangled fixed point of our RG transformation: the 
correlation matrices {r^^\r*^^\ ■ ■ ■ } quickly converge to a fixed r*;^^^. In particular, the 
renormalized entanglement of a block is constant. Thin lines correspond to a non-critical 
system, A = 1.001, which the RG fiow, as generated by ER, eventually brings a product 
(unentangled) ground state. Right: Quantum XX model, 7 = 0. Bold/thin lines represent 
two critical cases, A = and A = cos(157r/16). They belong to the same universality class 
and are found to indeed converge to the same correlation matrix r*^, (with r*^ 7^ ^ising) 
and in particular to the same renormalized entropy. 



is remarkably small: the error in the critical ground state energy is less than 10~^, while 
the two-point correlators (aja^), reconstructed from the MERA, accumulate a relative 
error that ranges from 10"'' for nearest neighbors to 10% for |r — s| ~ 4,000. Had we 
not used disentanglers, the error in the energy would be 10~^ after only a single RG 
transformation and an error of 10% in the two-point correlators is already achieved for 
|r — s| = 42. (ii) The line 7 = [equivalent to the quantum spin XX model] is critical 
for |A| < 1. Here we consider P = 4 modes per site and apply again 13 iterations of the 
RG transformation, reaching sizes of 4 x 2^^ = 32768 modes. The errors in energy and 
correlators are similar to those in the line 7 = 1. In both cases, an analysis of the RG 



flow and its fixed points in terms of entanglement is quite insightful, see Fig. |3.2[ ER can 
also be used to generate a RG transformation in the space of Hamiltonians, by replacing 



Eq. 3.12 with a minimization of the energy. Fig. 3.3 shows that critical systems are also 
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Figure 3.3: Dispersion relation of Hamiltonian (3.1) in ID with 7 = 0, the quantum 
spin XX model, under successive RG transformations. Shading indicates the Fermi sea. 
A sequence of local, coarse-grained Hamiltonians is obtained {H^^\ H^'^\ ■ ■ ■ } with their 
corresponding dispersion relations {1^1,1^2, '} converging to a straight line, a fixed point 
of the RG flow. Convergence is achieved very quickly at half filling (A = 0) and slower for 
A = cos(37r/4). These results have been obtained by minimizing the energy while keeping 
8 modes in each effective site. 



fixed points of this alternative approach, that preserves the low energy spectrum. 

In 2D the model has three phases, denoted I, II and III in Ref. |LDY"'"06] . where the 
distinct forms of entanglement scaling were characterized. In phases II (critical, with a 
Fermi surface consisting of a finite number of points) and III (non-critical, with a gap 
in the energy spectrum) we are once more able to coarse-grain the system in a quasi- 
exact, sustainable manner. This is remarkable. The entropy of a square block made 
of modes grows as the size of its boundary, Sl ^ L |LDY"'"06] . This implies that 
the number of modes we should keep in an effective site grows exponentially with the 
number of iterations of the RG transformation, which is precisely why DMRG does not 
work for large 2D systems. Instead, disentanglers bring this number again down to just 
a constant. As a result one can, in principle, explore systems of arbitrary sizes. In 
particular, by considering P = 4^ modes per site we apply r = 4 iterations of the RG 
transformation, with a final block effectively spanning P x 4^+^ = 16384 modes, whilst 
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Figure 3.4: Entanglement entropy Sl of a block oi L x L modes in 2D models. Left: 
in the critical phase II and the non-critical phase III (bold/fine lines respectively) the 
entanglement entropy grows linearly with the size L of the boundary of the block, Sl ^ L 
(boundary law). As in ID, the renormalized entanglement is constant for the critical 
model and it eventually vanishes for the non-critical model. We have considered 7 = 1 
and A = 2, A = 2.05 for the critical/non-critical case. Right: The critical phase II system 
(7, A) = (1, 2) is replotted for comparison against critical phase I, (7, A) = (0, 0), where the 
system has a ID Fermi surface and the entanglement entropy has a logarithmic correction, 
Sl ~ LlogL. Here disentanglers are not able to reduce the renormalized entanglement 
down to a constant. 



maintaining truncation errors of the same scale as the ID models analyzed, = 1.1 x 10~^. 
As in the ID case, the structure of fixed points of the RG flow can be understood in terms 
of the renormalized entanglement, see Fig. 3.4 On the other hand. Phase I (critical, with 
a one-dimensional Fermi surface) is so entangled that ER is no longer able to prevent the 
growth in the number P' of modes that need to be kept per site. The system displays 
a logarithmic correction to the entropy, Sl ~ LlogL |LDY+n6l IB(]Sn6[ IWolOBi [GK06], 
while the MERA can only reproduce a linear scaling Sl ^ L |Vid08j if just a constant 
number of modes are kept per site, P' = P. 
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3.4 Conclusions 



We have presented, in the simphfied context of fermion models with quadratic Hamil- 
tonian, unambiguous evidence of the vahdity of the ER approach in ID and 2D sys- 
tems. Similar derivations can be also conducted for bosonic lattice systems with quadratic 
Hamiltonians, as is investigated in Chap. |4| These results show that the MERA |Vid08j 
is an efficient description of certain 2D ground states. A number of examples also con- 
firm that [i) ER produces a quasi-exact, real-space RG transformation where the coarse- 
grained lattice is locally equivalent to the original one, enabling the study of RG fiow 
both in the space of ground states and Hamiltonians; (ii) non-critical systems end up in 
a stable fixed point of this RG fiow, where the corresponding ground state is a product 
(i.e. fully disentangled) state, whereas scale invariant critical systems end up in an un- 
stable fixed point, with an entangled ground state. Moreover, ER sheds new light into 
the ground state structure of two-dimensional systems with a one-dimensional Fermi sur- 
face. There, the presence of logarithmic corrections in the entropy Sl of a large block 
|LDY+06[ IBCS061 IW0IO6I IGK06] cannot be accounted for with the MERA with constant 
number of modes P per level of coarse-graining, hinting for the need for a generalized 
MERA in order to properly describe such systems. 



Chapter 4 

Entanglement renormalization in free 
bosonic systems: real-space versus 
momentum-space renormalization group 
transforms 



4.1 Introduction 

The renormalization group (RG) is a set of tools and ideas used to investigate how the 
physics of an extended system changes with the scale of observation [KGH+GTI IWilTSj 
[Fii98l ICar96l IShi99l iDelOil ISha94l IWhi92l IWhi93l ISchOSj . The RG plays a prominent 
role in the conceptual foundation of several areas of physics concerned with systems that 
consists of many interacting degrees of freedom, as is the case of quantum field theory, 
statistical mechanics and condensed matter theory [KGH+G?! IFis98t ICar96t IShi99l IDel04l 



ISha94j ■ In addition it also provides the basis for important numerical approaches to study 
such systems |Wil75| IWhi92| IWhi93l ISchOSj IMW941 IMW96j . 

Given a microscopic description of an extended system in terms of its basic degrees of 
freedom and their interactions, RG methods aim to obtain an effective theory, one that 
retains only some of these degrees of freedom but is nevertheless still able to reproduce 
its low energy (or long distance) physics. The effective theory is obtained through coarse- 
graining transformations that remove those degrees of freedom deemed to be frozen at the 
observation scale of interest. For instance, given a Hamiltonian for an extended sys- 
tem one may aim to use successive RG transforms to obtain a sequence of coarse-grained 
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Hamiltonians (^H^^\ H^^\ H'^'^\ . . each an effective theory describing the original sys- 
tem at successively lower energy scales, c.f. Fig. Broadly speaking, RG techniques 
fall into two categories depending on how the coarse-graining is implemented, namely 
momentum-space RG |Dein4] and real-space RG |Whi92[ IWhi93[ ISchOSi IMW941IMW96] . 

Momentum- space RG is applied to theories that are expressed in Fourier space. It works 
by integrating out high-momentum modes of a field and it is often associated to pertur- 
bative approaches. Instead, real-space RG is applied directly to theories that are written 
in terms of local degrees of freedom, say spins in the case of a spin system defined on 
a lattice. It is not linked to perturbation theory and can in particular be applied to 
strongly interacting systems. As proposed by Kadanov |KGH"'"67] . the coarse-graining 
transformation is implemented by replacing a block of spins with a single effective spin, a 
procedure refined by Wilson |Wil75j and subsequently turned by White into the density 
matrix renormalization group (DMRG) algorithm |Whi92t [Whi93t ISchOSj . an impressively 
precise numerical tool to study one-dimensional systems. 




Figure 4.1: Given a description of an extended system in terms of a Hamiltonian 
one may use successive RG transforms T to obtain a sequence of Hamiltonians 
(^H^^\ H^^\ H^'^\ . . ^ each an effective theory describing the original system at succes- 
sively lower energy scales (or longer distances). Characterising the fixed points of the RG 
flow may help provide an understanding e.g. of the low-energy behavior of the system in 
the thermodynamic limit or of the stability of the system under various perturbations. 
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A major difficulty of momentum-space RG comes precisely from the fact that it requires, 
as a starting point, a description of the system in Fourier space. Such description is not 
always available and might not be obtained easily. Consider for instance a system of 
interacting spins, as specified by some generic spin-spin interaction. There, obtaining a 
momentum-space representation might be as difficult as solving the whole theory. In this 
and many other RG approach must be performed in real space. 

In spite of its indisputable success, the DMRG algorithm |Whi92l IWhi93t ISchOSj suffers 
from a shortcoming that has important implications. Because of the accumulation of 
short-ranged entanglement near the boundary of a spin block, the dimension of the Hilbert 
space used to effectively describe the block must grow with each iteration of the RG 
transformation. As a result, for instance, unstable fixed points of the RG flow (scale 
invariant critical systems) cannot be fixed points of the DMRG algorithm. Another, more 
practical consequence of this growth is that it limits the size of ID critical systems that can 
be analyzed and, most importantly, it severely limits the success of DMRG computations 
in higher spatial dimensions. 

Entanglement renormalization (ER) is a real-space RG method proposed in order to over- 
come the above difficulties |Vid07] . The main feature of ER is the use of disentanglers. 
These are unitary transformations, locally applied near the boundary of a spin block, 
that remove short-ranged entanglement before the system is coarse-grained. As a result, 
the effective dimension of the Hilbert space for a spin block can be kept constant un- 
der successive RG transformations, so that the approach can be applied to arbitrarily 
large systems. In this Chapter we explore the ability of entanglement renormalization 
to produce a sensible RG flow, one with the expected structure of fixed points and flow 
directions according to momentum-space RG, in D = 1,2 dimensional harmonic lattice 
systems. Such systems are an ideal testing-ground for ER. On the one hand, they have 
well studied properties |PEDC05l IAEPW021 ISkrOSl ICEPiOGj and can be fully character- 
ized in terms of correlation matrices, a fact that simplifies the analysis and conveniently 
reduces the computational complexity of ER calculations. On the other hand, an RG 
analysis of free-particle theories can be conducted simply and without approximations in 
momentum-space allowing for a comparison between the numerical results obtained using 
ER and the exact solution. The setting of harmonic lattices also allows for ER to be 
formulated in the language of bosonic modes, a formalism familiar to researchers in the 
of condensed matter physics and quantum field theory. 
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The results of this Chapter, presented in Sect. |4.5[ demonstrate that a real-space RG 
transform based upon entanglement renormalization is able to reproduce the exact results 
from momentum-space RG to a high accuracy in D = 1,2 dimensional lattice systems, 
both for the critical and non-critical cases considered. Also demonstrated is the ability of 
the MERA to provide an efficient and accurate representation of the ground state of free 
boson systems, thus extending to the bosonic case the results of Ref. |EV10bj and those 
of the previous Chapter. These results provide strong evidence that the ER approach, 
which can be implemented without making use of the special properties of free-particle 
systems, could be used to investigate low-energy properties of strongly interacting systems 
not tractable with momentum-space RG approaches |DEO08[ IRMV08[ lEVOQaj . 



The Chapter is organized in sections as follows (see also Fig. 4.2). Sect. 4.2 introduces 



the harmonic systems under consideration. In Sect. 4^ the process of renormalizing 
the system in momentum-space is explained, highlighting conceptual features of the RG. 
The critical system and examples of relevant and irrelevant perturbations are considered. 



Sect. 4.4 explains the details of the real-space RG implementation, both in terms of 



renormalizing the Hamiltonian and in terms of renormalizing the ground state directly. 



In Sect. 4.5 a comparison between results obtained from exact momentum-space RG and 



the results from numerical real-space RG is presented. 



4.2 Coupled harmonic oscillators 



In this Chapter the low-energy subspaces of harmonic lattices in D = 1, 2 spatial dimen- 
sions are to be analysed using both real-space and momentum-space RG transformations. 
We begin with a brief introduction to harmonic lattices detailing equivalent represen- 
tations in real-space and momentum-space coordinates, and the Fourier transform that 
shifts between such descriptions. For clarity the following derivations shall only be pre- 
sented for the ID system as the generalization to 2D, or higher dimensional systems, is 
straight forward. The Hamiltonian for a chain of harmonic oscillators each with mass 
m, angular frequency u, and coupled with nearest neighbors via 'spring constant' K, is 
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Figure 4.2: An outline of this Chapter. The harmonic lattice Hamiltonian Hus-, de- 
fined in terms of interaction between local degrees of freedom, may be coarse-grained 
directly with a numeric implementation of a real-space RG transform, such as entangle- 
ment renormalization (ER), as outlined in Sect. 4.4 Iterating the RG transform r times 



we get the t^^ effective Hamiltonian -ffjj2. An effective Hamiltonian may also be ob- 



tained by first transforming the Hamiltonian to a momentum-space representation i^MS; 

Momentum- 



4.2 



4.3 



The 



via Fourier transform of the canonical coordinates, as described in Sect, 
space RG transforms may be applied (analytically) to i^MS as described Sect, 
dispersion relations from the real-space Hamiltonians H^^^ are compared t o th ose of the 
corresponding momentum-space Hamiltonians, 



if^g in the results of Sect. 



4.5 



written 

N 



H = Co2_^\ —p^ + -^Qr + K [qr+i - qr) 1 

r=l V / 

N 

= Y.[fr + mWql + 2k - g,)') (4-1) 

r=l 

where in the second line we have chosen cq = 2m and defined K = mK for convenience. 
Note that periodic boundary conditions are assumed. The operators pi and qi are the usual 
canonical coordinates with commutation [pkiqi] = i^Ski- In our present considerations it 
is convenient to focus on the critical (massless) Hamiltonian 

N 

Ho = Yl (Pr + 2^ - qrf) , (4.2) 



r=l 



though the non-zero mass case will later be reintroduced in Sect. 4.3.1 As a preliminary 
to the momentum-space RG, the Hamiltonian Hq shall be recast into momentum-space 
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Figure 4.3: A comparison of an RG iteration for a D = 1 dimensional system (left) 
in terms of a dispersion relation, E[k), in momentum-space and (right) in terms of a 
lattice in real-space. The RG transformation maps the r**^ effective theory H^'^'> into a 
new effective theory while preserving the low-energy (or long distance) physics of 

the original theory. The three steps which comprise the iteration are described in detail 
in Sect. 4.3| for momentum-space RG and Sect. 4.4. 1| for the real-space RG. There are 
many possible ways the real-space coarse-graining step (i) can be implemented, Fig. 4.4 
describes two non-equivalent implementations. 
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variables via Fourier transform of the canonical coordinates. The Fourier-space coordi- 
nates p and q are defined 



TV 

i ' 



N 

2TTirK/N 



r=l 

^Y^q^e-^^^'-^"' . (4.3) 



r=l 



Substitution of the Fourier-space coordinates brings the Hamiltonian of Eq. 4.2 into 
diagonal form, here a set of N uncoupled oscillators 



(7V-l)/2 

^0= E (p^. + 8irsin^(^)g^). (4.4) 

K=-(Ar-l)/2 



Defining k = 27rK,/{aN), with constant 'a' representative of the lattice spacing, the ther- 
modynamic limit {N — 7- oo) can be taken 

n/a 

Ho= j (^{kf + Sksm^[^^q{kf^dk. (4.5) 

fe=— 7r/a 

The theory has a natural momentum cut-off A = ±7r/a originating from the finite lattice 
spacing; taking the lattice spacing 'a' to zero recovers the continuum limit with corre- 
sponds to the field theory of a real, massless scalar field. Two equivalent representations 
of the harmonic chain have been obtained; that of Eq. |4.2 written in terms of spatial 



modes (amenable to numeric, real-space RG) and that of Eq. 4.5 written in terms of 
momentum modes (amenable to analytic, momentum-space RG). The dispersion relation 
Eolk) of the system, which describes the energy of momentum mode k and is known from 
the solution to a single oscillator, is given by 

Eo{k) =2V2K\sm{ka/2)\. (4.6) 



The RG transformations of the Hamiltonian shall be chosen such that the resulting ef- 
fective theory preserves the low-energy structure of the original theory, or equivalently, 
preserves the small k part of the dispersion relation. As the low-energy part of the dis- 
persion is gapless and linear we would expect that, under the RG flow, the renormalized 
Hamiltonians would tend to a fixed point that has a linear, gapless dispersion. 
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4.3 Momentum- space RG 



In this section the harmonic system of Eq. 4.5, which has been recast in a Fourier basis 



is analysed using momentum-space RG. Although the RG analysis of interacting systems 
is more complicated that the free particle analysis undertaken here, and often involves 
perturbation theory, the basic procedure is the same. It is useful to consider the RG 
transformation as occuring in three steps, (i) Firstly the momentum cut-off is reduced, 
A (—)■ A' = A/2, and modes greater than the cut-off are integrated out of the theory. As 
there is no interaction between the momentum modes of Eq. |4.5 this step is presently very 



simple; the cut-off is reduced A A' = A/2 whilst leaving the form of the Hamiltonian 
for modes with momentum k < A' unchanged 

7r/2a 

Ho= j (^{kf + Sksm^{^^q{kf^dk. (4.7) 

fe=— 7r/2a 

(ii) Next the length associated to the system is changed, this can be implemented as a 
scaling of the lattice spacing a \^ a' = 2a, which gives 

Tx/a' 

H'^= j (^{kf + {^^q{kf^dk. (4.8) 

fc=— Tr/a' 

A change has been made to the observation scale of the system in terms of length. Next a 
change in the observation scale is made in terms of energy. Indeed, in the final step (iii) 
the fields are rescaled 

p{k) ^p'{k) = -^PW 

q{k)^q'{k) = ^q{k), (4.9) 

so that the new field operators have a modified commutation relation [p'(fc), q'ik)] = ih/2, 
in accordance with the desired change of energy scale. In principle the RG transformation 
is complete, however in this instance a further transform is required to recast the critical 
Hamiltonian into a manifestly invariant form. The field operators are rescaled once more 

p' ^ p" = V2p' 

q'^q" = ^q'. (4.10) 



^The approach of rescaling the lattice spacing is common in the context of condensed matter problems; 
equivalently we could have rescaled the momentum of the theory, k ^ k' = 2k, as is the approach most 
often used for the RG in a quantum field theory setting. 
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In contrast to the previous transform of Eq. 4.9, this transform is commutation preserving 



hence does not affect the physics of the system; namely the dispersion relation remains 
unchanged. Implementing the third step, together with the auxiliary transform of Eq. 



4.10 



the first renormalized Hamiltonian ^q^'' is given 

7r/a' 

Hj,'^ = J (^p" {kf + 32irsin2 q" {kf^ dk. (4.11) 

fc=— 7r/a' 

The RG transformation is summarized: (i) the degrees of freedom that are not relevant to 
the low energy physics are removed, followed by a changes of observation scale in terms 
of (ii) length (a i— )■ a' = 2a) and (iii) energy {h^ h' = h/2). Starting from a theory with 
a natural length scale, the lattice spacing a, the RG transform is thus used to derive an 
effective theory with the new length scale a'. The scale factors chosen in steps (ii) and 
(iii) may depend on the implementation of the RG as well as the problem to which it is 
being applied. Iterating the RG transform (dropping the 'primes' from notation), the r*^ 
renormalized Hamiltonian is 

ir/a 

iyM= I ^p{kf + 2'^^' K Sin' i^^^q{k?yk, (4.12) 

fc=— vr/a 

with corresponding dispersion relation 

E^;^ (k) = 2"+i I sin {ka/2^+^) \ . (4.13) 

In the limit of infinitely many transforms, r — )■ cxd, we get the Hamiltonian Hq°°'^ at the 
fixed point of the RG flow 

7r/a 

= J [p{kf + k{ea')q{kf'^dk. (4.14) 

fc=— 7r/a 

The fixed point Hamiltonian has a purely linear, gapless dispersion 

E^^\k) = a\/2k\k\. (4.15) 

as anticipated. 



4.3.1 Relevant perturbation 



In the previous section the massless harmonic lattice of Eq. 4.5 , a critical system, was 
shown to be a (non-trivial) fixed point of the RG flow. We now consider the stability of 
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this Hamiltonian under the addition of perturbations. Perturbations to the critical theory 
can be classified as being relevant or irrelevant depending on whether the deviations from 
the fixed point, induced by the perturbations, grow or diminish under the RG fiow. In 
order to study relevent perturbations a mass term H^e\ is reintroduced to the critical 



system of Eq. |4.5t as shall be shown shortly this term grows under the RG fiow. Hence, 
it significantly modifies the low-energy physics from that of the unperturbed system. The 
mass term is diagonal in both real-space and momentum-space representations 

n/a 

4ei = $^i'= J Hkfdk. (4.16) 

k=—iT/a 

The perturbed Hamiltonian H^n = Hq + m?H^e\ is equal to the original Hamiltonian 



of Eq. AA_ describing coupled harmonic oscillators of mass m. Since this perturbation 
term does not reorder mode energies {E{k) is still an increasing function of |A;|) the same 
RG transformations may be performed on Hj-^i as on the unperturbed system Hq. This 
analysis shows that the r**^ renormalized perturbation term H^^^ grows exponentially with 
the RG iteration r 

22-4,1, (4.17) 

Thus even the addition of even a small mass m to the critical system leads to a large 
difference between the perturbed Hm and original Hq Hamiltonians after only a few RG 
iterations; that is to say, the perturbed system tends to a different fixed point of the RG 

(t) 

fiow. This is further evidenced by consideration of the dispersion relation Em of the 
perturbed system H^. 



E^^^ (k) = 2"ym2 + S^sin^ {ka/2^+^) 

a^K , ^ ^ ( a?K . 
= 2^m+- — A;2 + — — - . 4.18 

In contrast to the linear, gapless dispersion e'^^ of the massless theory described by Eq 



4.15 we now see a quadratic dependence of the energy on the momentum fc, together 



with an energy gap that grows exponentially with the number r of RG iterations. 



4.3.2 Irrelevant perturbation 



Perturbations that become smaller along the RG fiow are termed irrelevant perturbations^ 
as these do not affect the asymptotic, low-energy behavior of the system. In this section 
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we construct and analyze an example of such a perturbation. The irrelevant perturbation 
ffirrei is constructcd from neighboring and next-to-nearest neighboring quadratic couplings 



Hirrel = ^[-Qr + (Qt+I - Qrf + T {qr+2 + qrf j • 
r=l ^ / 



(4.19) 



In the Fourier basis ifirrei has representation 

7r/a 

H'mcl 4 



sm 



y 



(4.20) 



We consider the perturbed system f^Q, = Hq + aHi„e\ with the perturbation strength 
chosen a > 0. The t^^ renormahzed perturbation -ff/^rei t)e obtained through the 
same sequence of RG transforms as applied to the unperturbed system 

Tv/a 



H.. 



irrel 



ka 



2'^^'sinU-^\q{krdk 



k=—TT /a 



■n/a 



a^k^ + O 



'k' 



)2t+2 



q {kY dk, 



(4.21) 



fc=— vr/a 



The perturbation Hnrei is thus exponentially suppressed under the RG flow. Equivalently, 
the addition of this term to the critical system has an effect on the low-energy physics 
that diminishes with each RG iteration, as is seen directly from the dispersion E^^ of the 
perturbed system Ha 



eM(A;)=2^+i 



sm 



ka 



'2K + asm^ 



ka 
2^+1 



= aV 2K \k\ + O {2-'^^) , (4.22) 

which converges to the same linear, gapless fixed point as did the unperturbed system of 
Eq. [415 



4.4 Real- space RG 

In this section an implementation of real-space RG based upon coarse-graining transfor- 
mations of the lattice is described. Following the seminal works of Migdal, Kadanoff, and 
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Figure 4.4: The coarse-graining step of the real-space RG, in which a block of two sites 
(each composed of M bosonic modes) of the original lattice C is mapped to a single site of 
the coarse-grained lattice can be accomplished in different ways. (Bottom) The sim- 
plest method involves applying a series of local projections, realized by isometric tensors 
w, whose function is to select and retain relevant block degrees of freedom. This shall 
subsequently be referred to as the local projection (LP) method. (Top) Entanglement 
Renormalization (ER) differs from the LP coarse-graining by the inclusion unitary disen- 
tanglers m, enacted across block boundaries, before the truncation with isometrics. Here 
we follow the binary coarse grainging scheme introduced in Chap. [2] 



Wilson in real-space RG |KGH+67l IWil75j . the coarse-graining transformation maps a 
block of sites from the original lattice C into a single site of a coarser lattice Let us 
consider a ID lattice C of N sites, each site described by a vector space V. We divide C 
into blocks of two sites and, following Wilson, implement a coarse-graining transformation 
by means of an isometry w 



w 



V'^V®^ w^w = Iy' (4.23) 



where V®^ is the vector space of two sites, V is the vector space of a site in the coarser 
lattice C of A^' = N/2 sites and ly is the identity in V. This coarse-graining trans- 
formation, which we shall refer to as a local projection (LP) transformation, is depicted 



graphically in Fig. 4.4 From an initial Hamiltonian H defined on lattice C we can obtain 



an effective Hamiltonian H' on lattice C via the transformation 



H' = W^HW, W = w^^/^. 



(4.24) 
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The LP transformation has the property of preserving locality of operators; for instance 



if the original Hamiltonian was a sum of two-body interactions, H = YliLi then the 

effective Hamiltonian would remain a sum of two-body interactions, H' = YliLi ^'i i+i- 
alternative coarse-graining transformation, known as entanglement renormalization (ER), 
follows as a modification of the LP scheme. As with the LP scheme, we map two sites 
of C into a single effective site of via an isometry w, however in ER one first enacts 
unitary disentanglers u 

y;V®2H^V®2, u^u = uv) = hm (4.25) 



across the boundaries of adjacent blocks, as show Fig. 4.4 Thus the effective Hamiltonian 



H' , as given by a transformation with entanglement renormalization, is 

H' = (lytf/t) H (f/^) ^ jj = u®N/2_ ^4 26) 

Inclusion of the disentanglers in the RG scheme, although having the unfortunate effect of 
increasing the computational cost of numeric implementation, has profound implications 
regarding the ability of the method to produce a sensible real-space RG flow, as shall be 



demonstrated in the results of Sect. 14.51 

Starting from the original Hamiltonian (C^^\ H^^^) = {C,H) one could, by iteration of 
the RG map, generate a sequence of effective Hamiltonians defined on successively coarser 
lattices 

^ (£(1), if (D) 4 (£(2)^ iy(2)^ . . . ^ (4.27) 

with the transformation T representing either a LP or ER transformation. The LP trans- 
formation T^p is characterised by the corresponding isometry w^'^\ while the ER trans- 
formation T^^^'' is characterised jointly by an isometry and a disentangler, {w^'^\u'''^^). 



Equation 4.27 can be used to directly investigate how H changes under scale transforma- 
tions allowing e.g. one to characterise the stability of H under perturbations or investigate 
properties of the system in the thermodynamic limit directly. 

In order to make meaningful comparisons between effective Hamiltonians defined at dif- 
ferent length scales, say between H^'^'^ and ii'(^+^) defined on lattices C^'^'^ and 
respectively, it is required that the dimension of the Hilbert space of a site in £("^+1) re- 



mains the same as that of C^^\ or equivalently that V = V in Eq. |4.23| This ensures 
that, for instance, two-body operators defined on C^'^^ and exist in the same pa- 

rameter space, allowing a direct comparison between the coefficients that describe the 
operators. Keeping the dimension of effective lattice sites in lattice C^'^'' constant between 
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RG transforms is also desirable for computational purposes. Indeed, if the dimension of 
effective sites were to grow with each RG iteration then the computational cost of imple- 
menting the transformations would also grow, limiting the number of transforms which 
may be performed. An investigation of whether the LP and ER transformations can 
accurately coarse-grain Hamiltonians over repeated RG iterations while keeping a fixed 
local dimension is a focus of this Chapter. 

Note that, in general, the disentanglers and isometrics {u,w) that best preserve the low- 
energy space of a Hamiltonian H^^'^ under coarse-graining will depend on the specific 
Hamiltonian H^^^ itself; in practice optimisation techniques such as |DEU08[ IRMVOSj 
lEVOQaj ■ or the techniques presented in Chap, [sj are required to compute the tensors 
{u, w). There are however some systems where analytic expressions for these tensors may 
be obtained jWMj iKRVnO] . 



4.4.1 RG of the harmonic lattice 



We now turn our attention to the analysis of the harmonic lattice system with ER[^ 
Following the ideas discussed in the previous section, several possible algorithms |DEO08i 
IRMVnSi lEVOQaj could be applied to study the low -energy subspace of the harmonic 
systems directly. Application of these algorithms requires only that the Hamiltonian in 
question is composed of a sum of local interaction terms, as is the case with the systems 
we consider. However, as this study is focused on Hamiltonians containing only quadratic 
couplings, this property can be exploited in order to significantly reduce the cost of the 



numerical RG as well as simplify the analysis of the results. The Hamiltonian of Eq. 4.1 



describing a ID harmonic chain of modes, can be written in a concise quadratic form 



2N 



(4.28) 



by defining a quadrature vector R= {p,q) with 

P2 



P = 



<<1 = 



\PN J 



92 



(4.29) 



\qN J 



^The LP coarse-graining may be viewed as a simplification of ER in which the disentanglers u are set 
to identity transforms. Therefore only implementation of ER need be explicitly described. 
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where "H, henceforth referred to as the Hamiltonian matrix, is a 2N x 2N Hermitian ma- 
trix. The coarse-graining transformations shall be chosen such that the effective Hamilto- 
nians also only contain quadratic couplings, described by some new Hamiltonian matrix 
"H'. Thus the RG analysis may be performed in the space of Hamiltonian matrices "H, as 
opposed to the (much larger) space of full-fledged Hamiltonians H and a more efficient 
realization of ER is possible. Retaining the quadratic form of the Hamiltonian entails 
limiting the disentanglers u and isometrics w which comprise the ER map to cannonical 
transformations, namely transformations that preserve commutation relations. Consider 
a transformation of the Hamiltonian matrix by a 2N x 2N matrix S 

U^H' = S^nS. (4.30) 

In order for the transformation to be commutation preserving it is required that the 
transform be a symplectic matrix, S G Sp(2A^, M). A symplectic transform can be 
characterized as leaving the symplectic matrix S invariant under conjugation, S'^T.S = S. 



Given our convention of grouping the quadrature vectors in Eq. 4.29 the symplectic matrix 
takes the form 

In 
In 



(4.31) 



with /tv as the N x N identity. Additionally, the Hamiltonians under consideration take 
even simpler form than Eq. |4.28[ as there is no coupling between p and q quadrature 
degrees of freedom in the harmonic chain, the Hamiltonian may be expressed as 

H = fp + fUqq. (4.32) 

It is thus convenient to restrict symplectic transforms S to those which preserve the form 



of Eq. 4.32; we only consider transforms of the type S = V ®V , with V a special orthog- 
onal transformation, V G SO(A^). It can be easily checked that V ®V is a. symplectic 
transform. In fact, this is an element of the maximal compact subgroup of Sp(2A^, M). 



The p-quadrature part of the Hamiltonian in Eq. 4.32 remains trivial under these trans- 
formations, allowing us to focus on the g-quadrature part l-iq of the Hamiltonian matrix, 
which transforms as 

= V^UgV. (4.33) 

Let us group a number M of contiguous bosonic modes of the ID harmonic chain together; 
each group of M modes shall henceforth be referred to as a site of the original lattice 
C. The disentanglers u, which act on two sites (that is, on 2M modes), are chosen as 
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special orthogonal transforms u G S0(2M). Isometries w are realized as a composition of 
a special orthogonal transform followed by a projection, w = woWproj., with 

Wo e S0(2M), 

^proj . = (Om©/m). (4.34) 

The transformation of the entire lattice is achieved by first constructing the direct sum 
of the local operators 

N/2M N/2M 

W = ^ w, f/ = n. (4.35) 

i=l 1=1 

Given the disentangler and isometry {u,w), the Hamiltonian Tiq is coarse-grained into a 
new Hamiltonian Ti^ defined as 

Uq' = W^U^ (Hq) UW. (4.36) 



By the definition of the isometry in Eq. |4.34 it is ensured that if the initial lattice C has M 



bosonic modes per lattice sites then the coarser lattices C^'^'^ also have M modes per lattice 



site. As discussed earlier in Sect. |4.4[ keeping the number of degrees of freedom per lattice 
site constant between RG maps is necessary to allow meaningful comparison of operators 
at different length scales. Note that the number of modes per site M plays the role of a 
refinement parameter; choice of a larger M retains more parameters in the description of 
the effective theory, yielding more accurate results, at the cost of greater computational 
expense. The simplest choice of a one-to-one correspondence between bosonic modes and 
lattice sites, i.e. setting M = 1, does not give sufficiently accurate numerics, hence the 
need for grouping M > 1 modes into each lattice site. 



4.36 



retains 



It is only for a proper choice of disentangler and isometry {u, w) that of Eq. 
the low-energy subspace of the original T-Lq. This proper choice is found by optimisation 
over all possible {u,w). It is desired that the RG transform be optimized to project 
onto the minimum energy subspace of the original Hamiltonian; a matrix equation which 
describes the minimization can be written 



mm 



HH'q}) , (4.37) 



with the effective Hamiltonian T-L'^ as Eq. 4.36 The matrix Ti'^ describes a Hamiltonian 
that is translation invariant between blocks of 2M modes, hence the trace of "H^ (which 
describes the entire system) may be minimised by minimising the trace of an individual 
block of H'g. An optimisation method, based upon alternating updates for isometries w 



and disentanglers u, can be used to find suitable {u, w) that minimise Eq. 4.37 and best 
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Figure 4.5: An iteration of entanglement renormalization, broken into three steps, is 
depicted in terms of the direct sum structure of the Hamiltonian matrix 1-Lq. Dark shaded 
squares in 1-Lq represent couplings within the site of M modes, light shaded squares are the 
couplings between sites (at most next-to-nearest neighbor). Step(i), removing 'fast' modes 
is realized by transforming 1-Lq by conjugation with the disentanglers u and isometries w. 
Step(ii), rescaling the momentum, is achieved by removing the zero rows/columns from 
H'q. Step(iii), rescaling the fields, is achieved by directly scaling liq by a factor of 2. These 

three steps combined take the r*^ renormalized Hamiltonian matrix "Hg^^ to the (r + l)*'^ 
renormalized Hamiltonian matrix 1-L'^q^^\ 



preserve the low-energy space of the original Hamiltonian Hq. The optimisation method 
to be used here is similar to the general algorithm described in Sect. IV of Ref. jEVOQaj . 



Assuming that suitable (m, w) have been obtained, the full real-space RG transformation 



of the Hamiltonian may be achieved in three steps, as illustrated Fig. 4.5, analogous 



to the three steps of momentum-space procedure described Sect. 4.3 Firstly, (i) the 



Hamiltonian matrix Hq is transformed with direct sums of the disentanglers and isometries 



as Eq. 4.36 Recall from Eq. 4.34 that an isometry w, which acts upon a block of 2M 
modes, consists of a special orthogonal transform followed by a projection, w = WoWproy 
The projection has form Wproj = (Om © Im) with the trivial (zero) part Om describing 
the M modes to be removed from the system and the identity part Im describing the M 
modes retained in the effective description. In the next step, (ii) the zero rows/columns, 
those which were acted upon in the previous step by Om, are removed from Hamiltonian 

(0)" 

matrix to form a new matrix Hq . (iii) The final step of rescaling the fields is realized 



by scaling the Hamiltonian matrix by a factor, the same rescaling as Eqs. |4.9| and |4.10 
for momentum-space RG is realized by defining "H, 



(1) 



2nf'^", with "HJ^^ as the first 
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renormalized Hamiltonian. Iterating this procedure r times gives the t renormahzed 
Hamiltonian 'H['^\ In the thermodynamic limit, — )■ oo, the transformation can be 
iterated an arbitrary number of times to obtain a sequence of Hamiltonian matrices 




each describing a theory with quadratic interactions between at most next-nearest-neighbor 
sites, and each defined on an identical lattice C that has M bosonic modes per lattice 
site. 



4.4.2 Ground state RG 



In the previous section an implementation of real-space RG that could be used to coarse- 
grain harmonic lattice Hamiltonians was described. We now detail a similar procedure 
which allows ground- states of the harmonic lattices to be coarse-grained directly. Since 
we are dealing with systems of free-particles, the covariance matrix 7 gives a complete 
description of the (Gaussian) ground state {ipcs)- In the case of a Hamiltonian as in Eq. 



4.32, the convariance matrix is of the form 7 = 7p © 7g, where 7^ and 7^ are defined as 

= '^{'^Gs\PiPj I^GS) 

(7,)^. = 2(^Gs|gig,|V'Gs). (4.39) 

The derivation of analytic expressions for (7p, 7^) is a standard calculation and can be 
found e.g. Ref. |AEPW02] , but is also presented in Appendix [B| for completeness. Again, 
we choose the disentanglers u and isometrics w to be canonical transformations. We 
obtain a sequence of increasingly coarse-grained states each defined by covariance matrix 

with 7'^°'' = 7 as the original ground state. The coarse-graining transformations of the 
ground state shall be realized by symplectic transforms S acting in the space of the 
covariance matrix, 7 1— )■ 7' = S'^'jS. As there is no correlation between p and q quadrature 
coordinates, i.e. 7 = 7p © 7g, we may restrict consideration to symplectic transforms 
5* G Sp(2A^, M) that are only of the form S = A® {A^'^)^ with A a real, invertible N x N 
invertible matrix. The covariance matrix 7 transforms under conjugation by matrix S, 
which implies that 

%^% = A^ (7p) A 

7, ^ 7; = (7,) {A-^f, (4.41) 
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yeilding a new state 7' = 7p © 7^- The disentanglers u, which act on two contiguous sites 
each of M modes, are thus reahzed as real, invertible 2M x 2M matrices that transform 



the covariance matrix as per Eq. |4.41[ Isometrics w, which act on two contiguous sites 
within a block, are realised as w = WoWpj-oj with wq as a real, invertible 2M x 2M matrix 
and ifproj = (Om © Im) a projection onto M modes of the block. The direct-sum of the 
local operators (m, w) is constructed 

N/2L N/2L 

W^proj = ^/^proj, W^O = ^0, (4-42) 

i=l i=l 

N/2L 

f/ = u. (4.43) 

i=l 

in order to coarse-grain the entire lattice. Starting from the covariance matrix 7 = 7p©7g 
describing the ground state of the original system and, for any choice of {u,w), a new 
state 7' = 7p © 7^ can be obtained on a coarser lattice C through an ER transform 

7; = {Wpro,W,'U-') 7, {iU-TiW,-TWp.oi) . (4.44) 

It is only for proper choice of tensors {u, w) that the above transformation correctly 
preserves the properties of the original state and produces a meaningful coarse-grained 
state. We now address the issue of how the proper tensors {u, w) can be found. In the 
application of the RG to the system Hamiltonian, both in momentum-space and real-space 
formulation, the modes truncated at each iteration were chosen as high-energy modes in 
order to leave the low-energy structure of the original theory intact. For the ground state 
RG a different criteria is required to judge which modes should be truncated from the 
system. The proper truncation criteria in order to preserve the ground state properties, 
proposed by White as part of his DMRG algorithm |Whi92l IWhi93] . requires that the 
truncation of a block should be chosen to keep the support of the density matrix for the 
block. In the present formulation of bosonic modes, in which the representation of the 
state is given by a covariance matrix as opposed to a density matrix, this rule imposes 
that only modes in a block that can be identified as being in a product state with the rest 
of the system can be truncated and safely removed from the description of the state. Thus 
in comparison with Hamiltonian RG, which was optimised to truncate modes such that 
effective theory had minimal energy, here we optimise for {u, w) so that the truncated 
modes have minimal entanglement with the rest of the system. 
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The entanglement of a block with the rest of the system is known to be related to the 
symplectic eigenvalues of 7 for the block |PEDC05| IAEPW02] . Let the covariance matrix 
(7)|b = {lp)\i3 © {lq)\B describe the correlations within a block B of two M-mode sites. 
The 2M symplectic eigenvalues Aj of the block B are the eigenvalues of the matrix formed 
by taking the product of matrices (7p)|i3 and (7g)|B 

A, = Spect{(7p)|e(7,)|B}. (4.45) 

The Heisenberg position-momentum uncertainty relation, which here may be simplified as 
(p^) (^^) ^1/4, enforces that all symplectic eigenvalues are positive and have magnitude 
greater than unity, Aj > 1 for all i. The entanglement entropy Si of a mode 'i' with 
symplectic eigenvalue Ai is 



Si 



(4.46) 



with f{x) = — xlogx. It is seen that the entropy Si is zero when mode i is in a minimum 
uncertainty state, Aj = 1, and an also that Si is increasing function of Aj. It follows 
that, if a mode with eigenvalue Aj = 1 can be identified within a block B, then we are 
assured that the mode is in a product state with the rest of the system and may be safely 
truncated. Hence the tensors {u,w) should be optimized to minimize the eigenvalues Aj, 
and thus the entanglement entropy, of the modes to be projected out. The projection 



PVproj in Eq. 4.42 was defined to select the modes to be retained; we now construct the 
complimentary projector W^proj 

N/2M 

^proj = (/2Af - Wproj) (4.47) 
1=1 

which projects onto the space of the modes to be truncated. The part 7 of the covariance 
matrix 7 that is to be projected out during the RG iteration may be written 7 = 7p © 7g 
with 

% = {KoiW.'U-') 7, {iU-TiW,~TW,roi) . (4.48) 

Given that the modes projected out of block B are to have minimum entropy, the trans- 
forms {u, w) should be chosen to minimise the matrix equation 

min(tr{(%)|B(7p)|B}). (4.49) 
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MERA 




TTN 




Figure 4.6: (Top) Applying successive RG maps based upon ER, as depicted in Fig. 4.4 



to the ground state leads to an approximate representation of the ground state in terms 
of a set of disentanglers and isometries {u^'^\ w^'^'') connected in a MERA tensor network. 



Here we utilise the binary MERA scheme introduced in Chap. [2j (bottom) Implementing 
the RG based upon an LP coarse-graining gives an approximation to the state in terms of 
a tree-tensor network (TTN) |SDV06l ITEV09] , a different class of ansatz for many-body 
states on a lattice. 
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As with the energy minimization described by Eq. 4.37 for the Hamiltonian RG , this 



equation is optimised variationally to find good disentanglers u and isometries w. 

The apphcation of RG transformations to the Hamiltonian could be interpreted as pro- 
ducing effective theories for the low-energy subspace of the original system. Is there a 
similar interpretation for the coarse-graining of the ground state? To address this question 
we assume that the ER map has been iterated r successive times on the original ground 
state 7^°^ to get the coarse-grained state 7*^'^-'. Each of the r coarse-grainings is charac- 
terised by a disentangler u and isometry w, thus the set of these transforms (u^^\w^^^) for 
j = 1,2, ... ,T characterise the sequence of RG maps. If the modes truncated during each 
iteration were in an exact product state, equivalently the eigenvalues of every mode 'i' 
truncated was Aj = 1, then the sequence of transformations could be inverted as follows. 
Starting from 7'^'^^ truncated modes (which were in a product state) are replaced back, 
and each of the transforms of Eq. 4.44| is inverted 



as to recover the exact original state 7'-'^''. The coarse-graining of the ground state can be 
interpreted as storing information about the short range properties of the state into the 
tensors (m*-*-*, w^*-*) while preserving the long range information about the state; the set of 
tensors (u^'^\w^^^) together with state 7'-^'' thus serve as a representation of the original 
state 7*^°-'. The set 'j^'^^ and form the multi-scale entanglement renormalization 



ansatz (MERA) |Vid08j . a variational ansatz for states on the lattice, c.f. Fig 4.6[ [If 
the ground state RG is performed with an LP coarse-graining, as opposed to one based 
upon ER, a tree tensor network (TTN) |SDV06l ITEV09] approximation to the state is 
obtained in terms of the isometries w*-*^]. For the harmonic lattices we consider, as with 
most non-trivial models, the coarse-graining cannot be performed exactly and a MERA 
will be an approximate, rather than exact, representation of the true ground state. In 
the present case, irrespective of how the transforms {u, w) are chosen, the modes to be 
truncated will still be slightly entangled as manifested in their symplectic eigenvalues, 
which will fulfill Aj > 1. This entanglement will be ignored during the coarse-graining 
thus limiting the precision with which the original ground state 7*-'^^ may be recovered. 

In the present setting of free bosons, in which the exact ground state 7*^°^ is already known, 
we compute the MERA and TTN approximations to the ground state via coarse-graining 
the ground state covariance matrix. 



The point of this exercise is that it allows us to assess the accuracy with which a MERA 
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Figure 4.7: Sequences of dispersion relations of Hamiltonians and after r = 0,1,3 RG 
transforms, comparing the real-space ER results (bold) with the exact momentum-space 
results (dashed, offset by 0.2). (Top Series) The critical Hamiltonian Hq tends to a linear 
dispersion under the RG map with ER, in good agreement with momentum-space results. 
(Middle Series) The addition of even small mass, m = 0.2, to the critical system gives a 
marked difference in the dispersion relation after r = 3 RG transforms. (Bottom Series) 
The effect of an irrelevant perturbation quickly diminishes under RG flow; by r = 3 

(3) (3) 

iterations the original and the perturbed dispersions are virtually identical, Ea ^ Eq . 



or a TTN can represent the many-body state. It also provides the opportunity to study 
the RG flow of the ground state 'y^'^^ under scale transformations. However, had our goal 
been to investigate properties of the unknown ground state of a system H with the help 
of real-space RG, it would have been absurd to assume knowledge of the exact ground 

V'gs / to the ground state 



state IV'gs) from the beginning. In that case, an approximation 

may be found through e.g. variational minimization of energy {ipcs 
[EV09aj or through an alternative method |DEO08llRMV08] . 



H 



tpGS ) as per Ref . 
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LP coarse-graining 




1 3 5 7 9 1 3 5 7 9 

RG Iteration 



Figure 4.8: The entanglement entropy S" of a site in the t coarse-grained ground state 
7^"^) of the ID harmonic chain. (Left) In the critical (massless) regime the entropy of a 
site is infinite, as realized by the infinite constant Q; however the change in entropy along 
the RG flow can be computed via a limiting process (see Appendix [B| . The entropy of 
the state renormalized with an LP coarse-graining increases by a constant with each RG 
iteration, reproducing the logarithmic growth law, Sl = (1/3) log2 L + c, as expected from 
conformal field theory jVLRK03t ICC04j , whilst the entropy of the system renormalized 
with ER remains constant along the RG flow. Furthermore the sequence of renormalized 
ground-states {7*-^'*, 7^^\ . . .} rapidly converge to a fixed state 7* under the ER transforms. 
(Right) For several values of finite mass, the entropy of the states renormalized with the 
LP scheme saturate at a length scale governed by correlation length, whilst the theories 
renormalized with ER factorize into a product state (zero entropy) at the approximately 
same length. 



4.5 Results and discussion 



4.5.1 RG flow of Hamiltonians 



In Sect. 4.3 the low energy subspace of harmonic chains were analysed exactly with 
momentum-space RG. This analysis yielded sequences of dispersion relations describing 
the RG flow towards a critical fixed point as well as the RG flow resulting from adding rel- 



evant and irrelevant perturbation terms, as given by Eqs. 4.15, 4.18 and 4.22 respectively 



The same harmonic systems can also be analysed numerically by applying successive real- 
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Initial Dispersion IstRG 2nd RG 3rd RG 




3i/a 7i/a 7i/a 7z/a 



Momentum, k 

Figure 4.9: (Left) Sequences of dispersion relations (in non-rescaled energy units) for the 
(top) gapped and (bottom) critical ID harmonic systems after r = 0,1,2,3 real-space 
RG transforms, comparing the results from the local projection (LP) method, ^'lP' 
those from the entanglement renormalization (ER) coarse-graining, E^-^. The numeric 
dispersions produced by ER agree with the exact results obtained from momentum-space 
RG (dashed) as to be almost visually indistinguishable; though small errors are noticeable 
near the momentum cut-off, k = ir/a. The dispersion relations given by coarse-graining 
performed with the LP scheme, whilst reasonably accurate after the 1*^* iteration, rapidly 
diverge from the exact result with successive RG iterations. In all cases the numeric RG 
was performed keeping M = 4 modes per site. 



space RG transformations T, based upon ER, following the variational method described 
This produces a sequence of Hamiltonian matrices (T-L^^\'H^^\'H^'^\ . . .), 



in Sect. 



4.4.1 



each an effective theory describing successively lower energy subspace of the original sys- 
tem. By construction, 'H^'^^ only contains quadratic couplings between nearest and next- 
to-nearest neighboring sites. The dispersion relations of the effective Hamiltonians H^^^ 
are found by Fourier transform in a similar manner as presented for the original Hamil- 
tonian in Sect. 4.2 Figure 4.7 shows the comparison of the analytic (momentum-space) 
and numeric (real-space) dispersion relations. 



The dispersion relations produced from numeric real-space RG with ER are seen to ap- 
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Initial Dispersion Initial Spectrum IstRG 3rd RG 




1 G 1 G 1 G 

Momentum Mode Number 

Figure 4.10: Sequences of energy spectra (in non-rescaled energy units) for the (top) 
gapped and (bottom) near-critical 2D harmonic lattice system after r = 0, 1, 3 real-space 
RG transforms. The performance of the numeric real-space methods, local projection (LP) 
and entanglement renormalization (ER), are benchmarked against the exact solutions 
from momentum-space RG (dashed). The energy spectra are obtained by sampling the 
dispersion relation E{ki, ^2) on a finite grid of G points, with G chosen very large, and 
then ordering the values, {Ei < E2 < ■ ■ ■ < Eq}- The spectra of the systems renormalized 
with entanglement renormalization, E^-^, retain a high level of accuracy through all r = 3 
RG iterations, though do lose precision near the high energy cut-off of the effective theory. 
Energy spectra E^Jp of the effective theories obtained with an LP coarse-graining have 
diverged considerably from the exact result by r = 3 RG iterations. The numeric RG was 
performed keeping M = 9 modes per site. 



proximate the exact results to a high degree of accuracy for the three of the cases consid- 
ered: a critical system, and relevant and irrelevant perturbations on the critical system. 
Furthermore, the critical Hamiltonian H^q^ converges to a fixed point of the RG flow, 
'Hq^'' = "H*, after r > 3 RG transforms. This fixed point Hamiltonian matrix "H* is 
manifestly invariant under further transformations, T{'H*) = Ti*, to within small numer- 
ical errors. Consequently the isometrics w and disentanglers u, which comprise the ER 
transform T, also converged to fixed points w^'^^ = w* and m*-^-* = u* for r > 3. The 
Hamiltonian matrix T-La^ of the critical system with the addition of an irrelevant per- 
turbation, H^a^ = + aHl^l,, converges to the same fixed point as the unperturbed 



4.5 Results and discussion 



61 



Hamiltonian Tia = = T-L* for r > 3 transforms. 



4.5.2 RG flow of ground states 



On the other hand, in Sect. 4.4.2 we describe a real-space RG to coarse-grain the ground 



state of the ID harmonic chain. The coarse-graining transformations can be based upon 



either the LP or ER schemes of Fig. |4.4[ The real-space RG produces a sequence of in- 
creasingly coarse-grained ground-states (7^''^ 7^^\ 7*^^\ . . .) , where 7*^''^ is the ground state 
covariance matrix of a harmonic chain with mass m. Fig. |4.8| displays the entanglement 



entropy 5*, as defined in Eq. 4.46 of a site in the coarse-grained state 7^"^^ as a function 
of the RG iteration r. 

The entropy of the ground state 7q^'* of the critical chain, m = 0, when coarse-grained 
with the LP scheme, increases by a constant with each RG iteration. More precisely, if we 
recall that a site in lattice C^'^'' corresponds to a block of L = 2^ sites of the original lattice, 
this growth of entropy reproduces the expected logarithmic scaling, Sl = (1/3) logg L + c, 
for ID critical systems jVLRK03| ICC04j . This growth demonstrates that it is impossible 
for a critical ground state to be a fixed point of the LP scheme. 

The ground state of a system with finite mass m can also be analysed with the LP 
method. Initially, the ground state displays the same logarithmic growth of entropy as 
in the critical case, but it saturates approximately after r = log2 C iterations of the RG 
transformation, where ( is the correlation length. 

Turning to the ER scheme, if disentanglers are included in the coarse-graining step then 
the entropy per site of state 70^"* remains constant under RG transforms. This is made 
possible by the disentanglers, which remove short-range entanglement at each iteration. 
Moreover, for the critical case, the sequence of successively coarse-grained ground-states 
7q^^ explicitly converges to a fixed point, 7q^'' = 7* for r > 3. The non-critical ground 
states 7*-'^'' of a harmonic chain with finite mass m converged to a trivial fixed point, a 
product state 7**, under RG iteration. This occurs at the length scale of the correlation 
length (. The numeric results have been obtained by keeping M = 4 bosonic modes per 
lattice site. The MERA approximation to the ground state, obtained from successive 
coarse-grainings with ER, proves remarkably accurate. For the critical system, which is 
the hardest to analyse from a computational point of view, the MERA obtained from 
r = 9 RG iterations can reproduce the exact correlators of the ground state to error 
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Figure 4.11: The mean energy -Ers, as defined Eq. 4.51 of tlie effective tlieories obtained 
through real-space RG, based upon either an LP (sohd) or ER (dashed) coarse-graining, 
are compared against exact momentum-space resuhs ^ms in terms of the relative error 
AE = (-Ers — Ems) / Ems- Performing the real-space coarse-graining with a larger number 
of modes per site M, which allows for more parameters to be kept in the description 
of the effective theories, is shown to give more accurate results for the real-space RG. 
(left) For the ID critical harmonic chain the mean of the energy spectrum obtained from 
renormalizing with ER, while keeping M = 4 modes per site, remains approximately 1% 
greater than exact value throughout the RG iterations. The spectra obtained from the LP 
method, even when using significantly larger M, only gave at best 10% accuracy after the 
same number of RG iterations, (right) For the 2D near-critical harmonic lattice, ER gives 
accuracies after 5 iterations of no less than 4% and 2% for M = 4 and M = 9 respectively; 
this is in contrast to the LP method which is only accurate to within 60% for M = 4 and 
30% for M = 9. Barring small fluctuations, the accuracy of the coarse-graining with ER 
remains relatively stable with RG iteration. 



bounded by 1 x 10 



4.5.3 Local projection vs entanglement renormalization 



Figs. 4^ and 4.10 present numeric dispersion relations for the harmonic lattices in D = 1 
and D = 2 spatial dimensions respectively, and compare results produced by coarse- 
graining with LP to those obtained by coarse-graining with ER. The dispersion relations 
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for the effective theories produced from the LP method are reasonably accurate after a 
small number of RG iterations. However, they rapidly diverge from the exact results 
in subsequent iterations. On the other hand, coarse-graining with ER is shown to keep 
significantly better support of the low energy subspace and, most importantly, to maintain 
accuracy over repeated RG iterations. The numeric spectra obtained with ER for the 2D 



lattices of Fig. |4.10| displays a loss of accuracy towards the high-energy cut-off, which 
indicates difficulty in keeping a sharp cut-off numerically. However, this is of little concern 
as the primary interest lies in the low-energy physics and not the high-energy cut-off of 
the effective theory. 

The average mode energy E of an effective theory, defined in terms of its dispersion 
relation E{k), is 

t Eik)dk (4.51) 



2A 



A 



for a ID system and similar for 2D. The average mode energy is used to qualitatively 
analyze the accuracy of the effective theories obtained with real-space RG, as presented 



in Fig. 4.11 Keeping more modes-per-site M, hence more parameters in the descrip- 
tion of the effective theory, increases the accuracy with which the numeric RG may be 
performed. However, even with the choice of a very large M, very large the LP method 
shows significant increase in error along the RG flow for both the ID critical and 2D 
near-critical harmonic systems. This figure also confirms what was established visually 



in Figs. 4.9 and 4.10 that coarse-graining with ER not only produces a more accurate 
low-energy theory than with LP, but also maintains constant accuracy over successive RG 
maps. 



4.6 Conclusions 

Real-space RG techniques, primarily in the form of the DMRG algorithm |Whi92t IWhi93i 
ISchOSj . have proved an invaluable tool for the numeric analysis of low-energy properties 
of extended quantum systems. However, this class of real-space RG methods (including 
the LP coarse-graining described here) suffer from an important deficiency: namely they 



do not reproduce proper RG flows. As demonstrated by Fig. 4^, a ID critical system 
could not possibly be a fixed point of the LP map due to the growth of entropy along the 
RG flow, which is caused by the accumulation of short-range degrees of freedom. This 
deficiency in turn limits the size of ID critical and 2D lattices that can be analysed with 
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such methods. 

In this Chapter we have demonstrated, by comparison against exact results from momentum- 
space RG, that a coarse-graining based upon entanglement renormalization can reproduce 
proper RG flows. Demonstrations included the analysis of a critical Hamiltonian "Ho and 
its ground state 70, which were shown to rapidly converge to non-trivial fixed points of 
the RG flow, and 7q respectively, and analysis of several non-critical systems which 
converged to trivial fixed-points of the RG flow. By addition of an irrelevant perturbation 
to the critical Hamiltonian "Ho, a new Hamiltonian T-ia was obtained that described the 
same phase as "Ho, but differed in the local interaction. As is required of a proper RG 
flow, "Hq, and "Hq converged to the same fixed point "Hq under the RG map defined by 
ER. It is important to note that the tensors (u, w), which comprised the real-space RG 
transform, were not chosen based on heuristic arguments or the desire for a particular 
outcome. Instead they were found through optimisation based upon either energy mini- 



mization, as in Eq. 4.37 or attempting to retain the support of the ground state, as in 



Eq. 4.49 



In addition, the coarse-graining transformation based upon ER was shown to induce a 
sustainable RG map; one that could be applied arbitrarily many times without significant 
loss of accuracy and without the need to increase the local dimension M of the effective 
theories (so that the computational cost is also kept constant). The sustainability of 
the ER based RG map allows investigation of low-energy properties in arbitrarily large or 
infinite ID and 2D lattices, as has also been demonstrated in recent studies in which local 
observables are evaluated directly in the thermodynamic limit |EV09al lEVOQbt lEVlOcj 
and critical exponents computed |PEV09| IMRGF09] without the need for finite size scaling 
techniques. 

The implementation of entanglement renormalization presented in this Chapter exploited 
properties of free-particle systems in order to reduce the computational cost. More general 
algorithms exist, as shall be described in the following Chapter, which allow implemen- 
tation without making use of the special properties of free particle systems. Thus it is 
possible to use real-space RG, based upon ER, as a means of investigating the low-energy 
properties of strongly-correlated systems where perturbative approaches are not valid. 
However, the implementation of ER in the interacting case is computationally expensive- 
especially so for 2D lattices. The application of ER to interacting quantum systems on 
2D lattices is investigated in Chap. [7] and Chap. [8j 



Chapter 5 

Algorithms for entanglement 
renormalization 



5.1 Introduction 

Entanglement renormalization provides a conceptual framework for implementing real- 
space RG transforms on lattice systems. As was described Chapter [2| ER is naturally 
related to a class of quantum many-body states, the so-called multiscale entanglement 
renormalization ansatz (MERA). Perhaps the most important application of ER, and 
certainly the application that is the focus of this thesis, is as a numerical tool for investi- 
gation of low-energy subspaces of quantum many-body systems. However, the usefulness 
of ER extends beyond that of solely as a numerical tool |AV08t IKRVOQt ISwi09] . Chapters 
[3] and |4] achieved the first step in realizing ER as a numeric tool for quantum many-body 
systems by demonstrating the ability of the MERA to accurately represent a variety 
of ground state and low energy wavefunctions in = 1,2 dimensional systems of free 
fermions and free bosons. Nevertheless, the variational methods used in these Chapters 
to optimise the MERA was specific to the free-particle systems under consideration. In 
this Chapter we present opimisations algorithms that allow one to compute the MERA for 
the ground-state or low-energy subspace of a generic local Hamiltonian. We also describe 
how expected values may be computed from the MERA and provide benchmark results 
for a variety of ID spin models. 



The Chapter is organised as follows. Sect. 5.2 explains how to compute expected values 



of local observables and two-point correlators. Central to this discussion is the past causal 
cone of a small block of lattice sites and the ascending and descending superoperators, 
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which can be used to move local observables and density matrices up and down the causal 



cone. Sect. |5.3| considers how to optimize a single tensor of the MERA during an energy 
minimization. This optimization involves linearizing a quadratic cost function for the 
(isometric) tensor, and computing its environment. In Sect. 5.4 we describe algorithms 
to minimize the energy of the state/subspace represented by a MERA. A highlight of the 
algorithms is their computational cost. For an inhomogeneous lattice with sites, the 
cost scales as 0{N), whereas for translation invariant systems it drops to just O(logA^). 
Other variations of the algorithm allow us to address infinite systems, and scale invariant 
systems (e.g. quantum critical systems), at a cost independent of N. Sect. 5.5 presents 
benchmark calculations for different ID quantum lattice models, namely Ising, 3-level 
Potts, XX and Heisenberg models. We compute ground state energies, magnetizations 
and two-point correlators throughout the phase diagram of the Ising and Potts models, 
which includes a second order quantum phase transition. We find that, at the critical 
point of an infinite system, the error in the ground state energy decays exponentially 
with the refinement parameter x, whereas the two-point correlators remain accurate even 
at distances of millions of lattice sites. We then extract critical exponents from the 
order parameter and from two-point correlators. Finally, we also compute a MERA that 
includes the first excited state, from which the energy gap can be obtained and seen to 
vanish as at criticality. 



5.2 Computation of expected values of local observ- 
ables and correlators 



We begin this section by briefly recalling the MERA formalism presented in Sects. 2^, 2^ 
and 2.4 of Chap. |2} Let C denote a D-dimensional lattice made of N sites, where each site 
is described by a Hilbert space V of finite dimension d, so that = V®^. The MERA 
is an ansatz to describe certain pure states |\E') G of the lattice or, more generally, 
to describe a xt dimensional subspace Yu ^ V^. The MERA for the site lattice may 
be grouped into T ^ log N different layers, where each layer contains a row of isometrics 
w and a row of disentanglers u. We label these layers with an integer r = 1,2, ■■■T, 
with r = 1 for the lowest layer and with increasing values of r as we climb up the tensor 
network. The MERA implicitly defines a sequence of lattices 



Co^ Ci^ y Ct, 



(5.1) 
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where £o = is the original lattice, and where we can think of lattice as the result of 
coarse-graining lattice £r-i- 

Let 0^''''+^' denote a local observable defined on two contiguous sites [r, r+1] of L. In this 
section we explain how to compute the expected value 

(ot'^'^+^Vt/ = tr(o["'''+ilp). (5.2) 



Here P is a projector (see Eq. 2.7) onto the XT-dimensional subspace ¥[/ C represented 
by the MERA. For a rank xt = ^ MERA, representing a pure state G Yc, the above 
expression reduces to 

(o['-''^+i])^ = (^lol'^'^'+^l 1^). (5.3) 



Evaluating Eq. 5^ is necessary in order to extract physically relevant information from 
the MERA, such as e.g. the energy and magnetization in a spin system. In addition, the 
manipulations involved are also required as a central part of the optimization algorithms 
described in Sects. 15.31 and 15.41 The results of this section remain relevant even in cases 
where no optimization algorithm is required (for instance when an exact expression of the 
MERA is known jAVOSj IKRVOQj l 



As explained below, the expected value of Eq. 5.2 can be computed in a number of ways: 



By repeated use of the ascending superoperator A, the local operator 



mapped onto a coarse-grained operator ot on lattice Ct- Eq. |5.2| can then be 
evaluated as the trace of the coarse-grained operator ot, tr(o['''^'+^]p) = tr(oT). 

Alternatively, by repeated use of the descending superoperator V, a two-site reduced 



density matrix p['"'''+^] for lattice C is obtained. Eq. 5^ can then evaluated as 
tr(o[-'-+i]p) = tr(o['-''^+ilp["''-+i]). 

More generally, the ascending and descending superoperators A and I) can be used 
to compute an operator Ot and density matrix pr for the coarse-grained 



lattice Cr- Eq. 5.2 can then be evaluated as tT{o^'^'^~^^^P) = tr(oT 



[r',r'+l] [r',r'+l]^ 



First we introduce the ascending and descending superoperators A and V and explain in 



detail how to perform the computation of the expected value of Eq. |5.2[ Then we address 
also the computation of the expected value 

(0)v,=tr(OP), 0^5^o[^''-+i], (5.4) 
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where O is an operator that decomposes as a sum of local operators in £; as well as the 
computation of two-point correlators. Finally, we revisit these tasks in the presence of 
translation invariance and scale invariance. 

The ascending and descending superoperators are an essential part of the MERA formal- 
ism that was introduced in Ref. |Vid08j (see e.g. Fig. 4 of Ref. |Vid08j for an explicit 
representation of the descending superoperator V). These superoperators have been also 
called MERA quantum channel/MERA transfer matrix in Ref. |GMF08j . 



5.2.1 Ascending and descending superoperators 

In the previous section we have seen that the MERA defines a sequence of increasingly 
coarser lattices {Cq,Ci,-- - ,£t}- Under the coarse-graining transformation Ul of Eq. 
2.4, a local operator o^l'^^^\ supported on two consecutive sites [r, r-|-l] of lattice £^-1) is 
mapped onto another local operator Ot '"^ supported on two consecutive sites [r', r'+l] of 



lattice Cr (Fig- 2.3). This is so because in Ulo^^'I_i^^Ur most disentanglers and isometries 



of Ur and U]. are annihilated in pairs according to Eq. 2.2 The resulting transformation 



is implemented by means of the ascending superoperator A described in Fig. 5.1 

0[r',r'+l]^^(^[rv-+l])^ (5.5) 

In order to keep our notation simple, we do not specify on which lattice/sites the super- 
operator A is applied, even though A actually depends on r, r and r'. Instead, when 
necessary we will simply indicate which of its three structurally different forms (namely. 



left Al^ center Ac or right Ar in Fig. 5.1) is being used 



As above, let [r, r+1] denote two consecutive sites of lattice £r-i and let [r', r'+l] denote 
two consecutive sites of lattice Cr that lay inside the past causal cone of [r, r+1] G Cr-i- 



Given a density matrix pr """^^ in the descending superoperator V of Fig. 5.2 produces 
a density matrix p^^'^i^^ in £r-i, 

pl'T^=T>{pr^'^'^). (5.6) 

Notice that the descending superoperator V (which depends on r, r and r') is the dual 
of the ascending superoperator A, V = A* . Indeed, as can be checked in Fig. |5.3| by 
construction we have that, for any o^^^'^i^^ and pt 

tr (ol^ltMpY^''^'^)) = tr {Aiot':^'^)pr^^'^'^) . (5.7) 
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(a) {b) O. (c) O 




Ascending Superoperators: 



' J 



o 




± = 




C8 



Figure 5.1: The ascending superoperator A transforms a local operator o^_i of lattice 
Cr-i into a local operator Or of lattice £^ (for simplicity we omit the label [r,r + 1] that 
specifies the sites on which Or-i and Or are supported). Depending on the relative position 
between the support of Or-i and the closest disentangler, the operator can be lifted to lat- 
tice C-r in three different ways, indicated in the figure as (a), (b) and (c). Correspondingly, 
there are three structurally different forms of the ascending superoperator A, namely left 
Al, center Ac and right Ar, indicated as (a'), (b') and (c'). Notice that the figure com- 
pletely specifies the tensor network representation of the superoperator, which is written 
in terms of the relevant disentanglers and isometries (and their Hermitian conjugates). 
An explicit form for the average ascending superoperator A of Eq. 5.28 is obtained by 
averaging the above three tensor networks. 
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Pt-1 



Pt-1 



4 £ , 



Pr-l 



Descending Superoperator: 

(a') /Vi = V^ip,) 




k I 



[b')p^_, = Vc(pd (C) p,., = V^{p,) 





Figure 5.2: The descending superoperator V transforms a local density matrix of 
lattice Cr into a local density matrix p^-i of lattice £r-i- Depending on the relative 
position between the support of pr-i and the closest disentangler, the density matrix pr 
canbe lowered to lattice £r-i in three different ways, indicated in the figure as (a), (b) 
and (c). Correspondingly, there are three structurally different forms of the descending 
superoperator V, namely left V^, center Vc and right Pr, indicated as (a'), (b') and (c'). 
An explicit form for the average descending superoperator T) of Eq. 5.34 is obtained by 
averaging the above three tensor networks. 
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Correspondingly, there are also three structurally different forms of the descending super- 
operators, namely left T>l, center Vc and right in Fig. 5.2 



Duality: 





Figure 5.3: The ascending and descending superoperators, A and are dual to each 
other, see Eq. |5.7 This becomes evident by inspecting the above figure, where the 
superoperators are explicitly decomposed in terms of disentanglers and isometries. 



5.2.2 Evaluation of a two-site operator 



We can now proceed to compute the expected value tr(o['''''+^]p) of Eq. 5.2 from the 



MERA. This computation corresponds to contracting the tensor network depicted in the 



upper half of Fig. 5.4 



In a key first step, the contraction of the tensor network for tr(of^'''~^^lP) is significantly 



simplified by the fact that, by virtue of Eq. 2^, each isometric tensor outside the past 
causal cone of sites [r, r+1] G £ is annihilated by its Hermitian conjugate. As a result, we 
are left with a new tensor network that contains only (two copies of) the tensors in the 



causal cone, as represented in the second half of Fig. 5.4 Because the past causal cones 
in the MERA have a bounded width, this tensor network can now be contracted with 
a computational effort that grows with N just as O(logA^). One can proceed in several 
ways: 

Bottom-top. — In the bottom-top approach, we would start by contracting the indices 
of 0''''''+^] and the disentanglers and isometries of the first layer (r = 1) of the causal 
cone; then we would contract the indices of disentanglers and isometries of the second 
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Expected value of a local observable: 





Fig ure 5.4: (i) Tensor network corresponding to the expected value tr(o['''''"''^^P) of Eq. 
5.2 The two-site operator represented by a four-legged rectangle in the middle of 

the tensor network. The shaded region represents the past causal cone of sites r, r + 1 G C 
(ii) All isometric tensors that lay outside the past causal cone of sites r, r+1 G L annihilate 
and we are left with a simpler tensor network. 
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Top-tensor contractions: 




Figure 5.5: (i) The top tensor transforms a two-site operator ot-i defined on lattice 
Ct-1 into a one-site operator (a xt x Xt matrix) ot on the top of the MERA. (ii) The 
two-site density matrix px-i on lattice Ct-i is obtained through contraction of the top 
tensor with its conjugate. Notice that pr-i, as well as all pr, are normalized to have trace 
tr(Pr) = XT- 



layer (r = 2); and so on (Fig. 5.6). However, this corresponds to repeatedly applying 
the ascending superoperator A on Oq^''^'''^' = o^'''^"''^'. Therefore this is precisely how we 
proceed, obtaining a sequence of increasingly coarse-grained operators 

supported on lattices Cq, Ci, £2, ■ ■ ■ , and Ct respectively. Here, the xt x Xt matrix ot 



at the top of the MERA is obtained according to Fig. |5.5| and the expected value of Eq. 
5.2| corresponds to its trace, 



tr(o^r,r+iip^ = tr(oT). (5.9) 

Top-bottom. — In the top-bottom approach, we would instead start by contracting the 
indices of the tensors in the top layer (r = T) of the causal cone; then we would contract 



the indices of the tensors in the layer right below (r = T — 1); and so on (Fig. 5.7). 



However, that corresponds to first computing a density matrix pr-i for the two sites of 



Ct-1 according to Fig. |5.5| and then repeatedly applying the descending superoperator 
P. Therefore this is how we proceed, producing a sequence of two-site density matrices 

Pt-1 ■■■ P2 Pi Po [b.iUj 

supported on lattices £r-i, ■ ■ ■ , ^,2, Ci and Co respectively!^ The last density matrix 
p[r,r+i] ^ pIj^''""*"-'^] describes the state of the two sites of C on which the local operator o[''"'^'+^] 
is supported. Therefore we can evaluate the expected value of o^'^''^+'^\ 

tr{o^r,r+i]p^ = tr(o[^'''+ilp''''^'^)- (5.11) 



"'^Each operator Pt- in Eq. 5.10 is both Hcrmitian {p^. — p^-) and non-negative (((/)| Pr > 0, V0) but 
its trace is tripr) = Xt- For simplicity, we call pr a density matrix for any xt > Ij even though it is only 
a proper density matrix for xt — 1- 
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Bottom-top contraction: 
(i) 




Figure 5.6: The contraction of the tensor network in the lower half of Fig. 5A_ using the 
bottom-top approach corresponds to employing the ascending super operator A a number 
of times. In this particular case, we first use (i) Ar, then (ii) Ac and then (iii) Al, to 
bring the tensor network into a simple form whose contraction gives a complex number: 



the expected value of Eq. 5.2 
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Top-bottom contraction: 




Figure 5.7: The contraction of the tensor network in the lower half of Fig. 5.4 using 
the top-bottom approach corresponds to first implementing (i) a top tensor contraction 
followed by repeated application of the descending super operator P. Specifically, here 
we first use (ii) V^, then (iii) Vc and then (iv) P^, in order to compute the appropriate 
density matrix ^['"'''+^1 for two sites [r, r + 1] G £. With the density matrix pl'^'^'+^l we can 
finally compute (v) the expectation value tr(of'"''"^^lP) = tr(o^^''^^^^ p^^'''^'^^^) . 
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Middle ground. — More generally, one can also evaluate the expected value of Eq. 5.2 



through a mixed strategy where the ascending and descending superoperators are used 
to compute 
which fulfill 



to compute the operator ot^'"^^^^^ and density matrix p^T^'''^'*'^' supported on lattice 



tr(o[-''-+i]p) = tr(oi'^-^'-+^lpi'^''^+'^)- (5.12) 

In all the cases above, one needs to use the ascending/descending superoperators about 
T ^ logiV times, at a cost 0(x*), so that the total computational cost is O(x^logA^). 

5.2.3 Evaluation of a sum of two-site operators 

In order to compute the expected value 

(0)v, = tr(OP), O^Yl ^'"''^'^ (5-13) 

r 

of an operator O on £ that decomposes as the sum of two-site operators, we can write 

tr(OP) = J2 tr(o['^'"+^lp) (5.14) 

r 

and individually evaluate each contribution tr(ot''''"+^lP) by using e.g. the bottom-top 
strategy of the previous subsection, with a cost 0(x^A^log A^). However, by properly 
organizing the calculation, the cost of computing tr (OP) can be reduced to 0{x^N). We 
next describe how this is achieved. The strategy is closely related to the computation of 
expected values in the presence of translation invariance, as discussed later in this section. 
Again, there are several possible approaches: 

Bottom-top. — We consider the sequence of operators 

Oo ^ Oi ^ O2 ^ •■■ Ot, Oo = 0, (5.15) 
where the operator Or is the sum of local operators, 

0^= ^olr''^+^l. (5.16) 

Or-1 is obtained from Or-i by coarse-graining. Or = UlOr-iUr- Each local operator 
jg ^YiQ sum of three local operators from Or-i (see (a),(b) and (c) in Fig. 



5.1), which are lifted to Cr by the three different forms of the ascending superoperator. 
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Al, Ac and Ar. Since Or-i has local operators, is obtained from Or-i by 

using the ascending superoperator A only N/3'^~^ times. Then, since Y1t=o'^^^ ^ 2' ^^^^ 



means that the entire sequence of Eq. 5.15 requires using A only O(A^) times. Once Ot 
is obtained, the expected value of O follows from 



tr(OP) = tr(OT) 



(5.17) 



Top-bottom. — Here we consider instead the sequence of ensembles of density matrices 

Et-1 ■■■ E2 ^ El ^ Eo, (5.18) 

where E^- is an ensemble of the N/3'^ two-site density matrices pr'^^^^ supported on nearest 
neighbor sites of Cr, 

i?.^{p[^'^pr-^V--,pf/^^'^^} (5.19) 

From each density matrix in the ensemble E^- we can generate three density matrices in the 
ensemble -Er- 1 by applying the three different forms of the descending superoperator, V^, 



T>c and (see (a),(b) and (c) in Fig. 5.2). All the N/'y^^ density matrices of ensemble 

Et-_i can be obtained from density matrices of E^. in this way. Since X]t=o ^ ^' 

that by using the descending superoperator V only O(A^) times, we are able to compute 



all the density matrices in the sequence of ensembles of Eq. 5.18 Once the ensemble 
Eq = E has been obtained, 

E={p^^p[2'^---,p[^'^l}, (5.20) 
the expected value of O follows from 

tr(OP) = ^tr(o['-''-+iV^'''^'')- (5.21) 



Middle ground. — More generally, we could build operator Or as well as ensemble inl- 
and evaluate the expected value of O from the equality 

tr(OP) = 2];tr(o[r'^'+^Vi'''^'')- (5.22) 



Each of the strategies above require the use of the ascending/descending superoperators 
O(A^) times and therefore can indeed be accomplished with cost 0{x^N). 
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Figure 5.8: In order to compute a two-point correlator C2(ri, we need to consider the 
union of the past causal cones of sites ri and r2. Notice that, in contrast with the case of 
a single local operator, the joint causal cone of two distant sites typically involves more 
than two contiguous sites of some lattice This makes the computational cost scale as 
a power of x larger than x^. 



5.2.4 Evaluation of two-point correlators 

Let us now consider the computation of a two-point correlator of the form 

C2(ri, r2) = (^1 ol"^] ® o["^] 1^) , (5.23) 
where ot^' and o^^'^ denote two one-site operators applied on two arbitrary sites r and s 



of £, see Fig. 5.8, Fig. 5.9 shows the tensor network to be contracted. Again, we can 
use Eq. 2^ to remove all disentanglers and isometries that lay outside the joint past 
causal cone for sites r and s. Then, we can proceed to contract the resulting tensor 
network, for instance through a bottom-top or top-bottom approach, with the help of the 
ascending and descending superoperators (and generalizations thereof). Notice that since 
at intermediate layers the two legs of the causal cone may contain two sites each one, in 
general we will need to compute operators/density matrices that span more than just two 
sites, and the cost of their computation will be larger than 0(x^)- 

However, for specific choices of sites r, s G £, we are still able to compute C2(r, s) with 
overall cost O(x^logA^), as illustrated in Fig. 5.10 We emphasize that this was not 



possible in the binary ID MERA and is one of the main reasons to work with the ternary 
ID MERA. For such choices of sites r and s, each of the two legs of the joint past causal 
cone contains just one site until, at some layer tq, they fuse into a single two-site leg. We 



can introduce one-site ascending and descending superoperators A^^^ and 'D'-^-' (Fig. 5.11 ) 
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Two-point correlator: 




Figure 5.9: (i) Tensor network to be contracted in order to evaluate a two-point correlator 
C2{ri,r2). Similarly to the case of a local observable Fig. 2.3, the tensors outside of the 

The resulting 



2.2 



casual cone annihilate in pairs due to their isometric character, Eq 
tensor network (ii) is much simpler network. However, for a generic pair of sites ri, r2 G C, 
the joint past causal cone will contain more than just two sites per layer, resulting in a 
computational cost that scales with x as a power larger than x^- 
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in terms of which we can express, for r < Tq, the transformation of a product operator 
^ into a product operator 



or]®ori=^«(o!rii)®^«(oWi), 



or of a density matrix pi*" ''^ ^ into a density matrix 



(5.24) 



(5.25) 



where r, s G £r-i and r', s' G £r are sites corresponding to single- site legs of the causal 
cone. In, say, the bottom-top approach we can compute the correlator of Eq. |5.23| by 
using the single-site ascending superoperator A^^^ for layers r < tq and then the two-site 
ascending super-operator A for layer r > tq. 



(i) 




(iii) 
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Figure 5.10: Two-point correlators for specific pairs of sites r, s [at distances of 3'^ sites 
for g = 1, 2, 3...] can be computed with cost O(x^logiV). This is due to the fact that the 
causal cones for each of r, s contains only one site until they meet — (i) at £i, (ii) at £2 
or (iii) at £3. 



5.2.5 Translation invar iance 



The computation of the expected value tr(o[^'''"''^'P) of a single local operator 

the case of a translation invariant MERA (see Fig. 2.8) can proceed as explained earlier 
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(ii) One-site Ascending 
Superoperator: 



(iii) One-site Descending 
Superoperator: 




Pr-1 




Figure 5.11: A one-site operator o^-i supported on certain sites of Cr-i (corresponding 
to the central wire of an isometry Wr) is mapped onto a single-site operator on Cr- In 
this case the (i) ascending and (ii) descending superoperators A^^^ and T)^^^ have a very 
simple form. 



in this section. In the present case one would expect the result to be independent of the 
sites [r, r + 1] G C on which the operator is supported, but a finite bond dimension x 
typically introduces small space inhomogeneities in the reduced density matrix p['''''+^l and 
therefore also in tr(o['^'^+ilP) = tr(o[^'''+i]p^''''^^^')- 

Given a two-site operator o, an expected value that is independent of [r, r + 1] can be 
obtained by computing an average over sites, 

tr(oKr+i]p) ^ l^tr(o["''^+ilp) (5.26) 



i^5^tr(o[^''^+ilp'^'''+'^), (5.27) 



where the terms o^^'^'^^^ denote translations of the same operator o. This average can be 
computed e.g. by obtaining the N density matrices p['"''"+^l individually and then adding 
them together, with an overall cost 0{x^N). However, with a better organization of the 
calculation the cost can be reduced to 0(x*logA^). 



We first need to introduce average versions of the ascending and descending superoper- 
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ators. Given a two-site operator Or-i in lattice Cr-i, we can build a two-site operator 
Or by using an average of the three two-site operators resulting from lifting Or-i to lat- 
tice Cr, namely AL{or-i), Aiior-i) and ALior-i). In terms of the average ascending 
superoperator A, 

A=^iAL + Ac + AR), (5.28) 

this transformation reads 

Or = A{Or-l). (5.29) 

Importantly, if we coarse-grain the translation invariant operator 

^J^ot^l iV..^iV/3^ (5.30) 

r 

where A^^-i is the number of sites of £7-1 and the terms denote translations of Or-i, 

the resulting operator can be written as 

(5-31) 

^ r 

where the terms o^^^^^^ denote translations of Or and where Or and o^_i are related through 



Eq. 5.29 In other words, the average ascending superoperator A can also be used to 



characterize the coarse-graining, in the translation invariant case, of operators of the 



form of Eq. 5.13 



Let pr denote the two-site density matrix obtained by averaging over all density matrices 
p[r,r+i] different pairs [r, r + 1] of two contiguous sites of Cr, 

r 

and similarly for lattice Cr-i, 

(5-33) 

I\r—1 

r 

Recall that each density matrix pr'^'^^^ on lattice Cr gives rise to three density matrices 
in Cr^i according to the three versions of the descending superoperator, namely Vl, Vc 
and Vr. It follows that the density matrix pr^i can be obtained from the density matrix 
Pr by using the average descending superoperator, 

V = ^{Vl + Vc + Vr), (5.34) 
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that is 

Pr-l = A{pr). (5.35) 

We can now proceed to compute the average expected value of Eqs. 5.26 5.27 This can 
be accomplished in several alternative ways. 

Bottom-top. — Given a two-site operator o, we compute a sequence of increasingly 
coarse-grained operators 

oo A- oi A- 02 A- ■ ■ ■ ot, oq = o, (5.36) 

where Or is obtained from Or~i by means of the average ascending superoperator A. Then 
we simply have 

lj]tr(o['-'^+ilp)=tr(o^). (5.37) 

r 

Top-bottom. — Alternatively, we can compute the sequence of average density matrices 

Pt-1 ■■■ P2 -> Pi -> po, (5.38) 

where pr-i is obtained from by means of the average descending superoperator V and 
where p = Po corresponds to the average density matrix on lattice C, 

r 

in terms of which we can express the average expected value as 

l5^tr(o['-'^+V)=tr(op). (5.40) 



Middle ground. — As costumary, we can also use both A and V to compute Or and pr, 
and evaluate the average expected value as 

i5^tr(o[^''-+i]p)=tr(o.p.). (5.41) 

r 

In all the above strategies the average ascending and descending superoperators A and V 
are used 0(log(A^)) times and therefore the computational cost scales as 0(x*logA^). 



To summarize, with a translation invariant MERA we can coarse-grain a single two- 
site operator o (with a transformation that involves averaging over all possible causal 
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cones) or compute the average density matrix p by using the average ascending/ descending 
superoperators. This leads to a sequence of operators Or and density matrices pr, 

AAA 
Oo — )■ Oi — )■ ■ ■ ■ — )■ Ot, Oq = o, (5.42) 

Po ^ pi ^ ■ ■ ■ ^ Pt, Po = P, (5.43) 

from which the expected value of o is obtained as tr (op), as tr (ot) or, more generally, as 
tr (o^p^). 



5.2.6 Scale invar iance 



In the case of a translation invariant MERA that is also scale invariant (see Fig. 2.8), the 
average ascending superoperators A is identical on each layer r, since it is always made of 
the same disentangler u and isometry w. We then refer to it as the scaling superoperator 
S |PEV09j . Its dual S* corresponds to the descending superoperator V. 

As derived in Ref. |GMF08j . the expected value of a local observable o in the thermo- 
dynamic limit can be obtained by analyzing the spectral decomposition of the scaling 
superoperator S, 

S{») = ^ Aa0„tr(0a«), tr(0„0^) = 5ai3- (5.44) 

a 

We refer to the eigenoperators (pa of 5, 

S{(j)a) = KK, (5.45) 

as the scaling operators. Notice that the operators 0^, which are bi-orthonormal to the 
operators 0q,, are eigenoperators of iS*, 

S*{k) = Kk- (5.46) 



We recall that the scaling operator S is made of isometric tensors (cf. Eq. 2.2) and 
therefore the identity operator I is an eigenoperator of S with eigenvalue 1 (that is, S is 
unital), 

5(1) = I. (5.47) 

On the other hand, since the MERA is built as a quantum circuit — and descending 
through the causal cone corresponds to advancing in the time of a quantum evolution — 
it is obvious that the descending superoperator P is a quantum channel, and so are V 
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and S* (see also |GMF08j ). In particular, S* is a contractive superoperator |BR79j . which 



means that the eigenvalues Aq in Eq. 5.44 are constrained to fulfill |Aq,| < 1. In practical 



simulations |PEVn9j one finds that the identity operator I is the only eigenoperator of S 
with eigenvalue one, Aj = 1, and that | Ad < 1 for a 7^ I. Let p denote the corresponding 
unique fixed point of S*, 

S*{p) = p. (5.48) 

In an infinite system, the local density matrix of any lattice Cr (with finite r) results from 
applying S* on p^ an infinite number of times, and it is therefore equal to the fixed point 



p |GMF08j . Consequently, Eqs. 5.42 and 5.43 are then replaced with 

Oo = o, 



s s s 

. s* ^ s* ^ s* ^ 

p ^ p ^ p ^ p 



where in addition, by decomposing o in terms of the scaling operators (pa, 

O = ^ Ca(j)a, Ca = tr(0„o), 

a 

we can explicitly compute 0^. 

Or = ( t? O • • • O ^J (o) = ^C„(A„)^0„. 



(5.49) 
(5.50) 

(5.51) 
(5.52) 



This expression shows that, unless cj 7^ 0, the operator Or decreases exponentially with r 
(recall that |Aa| < 1 for a 7^ I) and its expected value must vanish. The average expected 
value of o then reads: 

l5^tr(o['-''-+i]p)=tr(op). (5.53) 



lim 



Two-point correlators for selected positions can also be expressed in a simple way, by con- 
sidering the one-site scaling superoperator S^^\ which is how we refer to the superoperator 



^(1) of Fig. 5.11 in the case of a scale invariant MERA. Its spectral decomposition, 

= ^(J'a'4'a^^i'i'a»), tT{^ai'p) = Sa)3, (5.54) 

a 

provides us with a new set of (one-site) scaling operators ipa- Given two arbitrary one-site 
operators o and o', we can always decompose them in terms of these ipa (similarly as in 



Eq. 5.51). Thus we can focus directly on a correlator of the form {ipa'^ilj^a^). Here r and s 



are restricted to selected positions as in Fig. 5.10 Then we have 
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where Aq, = logg/i^ is the scaling dimension of the scahng operator ipa, whereas Cap is 
given by 



Cai3 = tr((V^„(g)V^^)p), 



(5.56) 



with p from Eq. 5.48 



In deriving Eq. 5.55 we have used that, by construction, |r — s| = 3^ for some q 



1,2,3, 



. Coarse-graining ipa 'ip]^ a number q of times produces a muhiphcative factor 
{Ha^-iiY and the residual two-site operator ifj'ail'f^-, whose expected value gives Cap- On 
the other hand, by noting that /i'' = p'°S3k-*l = \r — s|^°§3M _ |^ _ ^jioggA^ arrive at 
{Hal^'pY = \r — sl^^^^f^ , which explains the denominator in Eq. 



5.55 



In conclusion, from the scale invariant MERA we can characterize the expected value 
of local observables and two-point correlators, as expressed in Eqs. |5.53| and 5.55 5.56 



All critical exponents of the theory can be extracted from the scaling dimensions A^. 
The scale invariant MERA is further explored in Chapter [6| including a description of 
an algorithm to optimise for a scale invariant MERA and benchmark computations of 
scaling dimensions Aq for several different spin models. 

The required manipulations include computing p from S (using sparse diagonalization 
techniques) and diagonalizing S^'^\ all of which can be accomplished with the ternary ID 
MERA with cost 0(x*). 



5.3 Optimization of a disentangler/isometry 

In preparation for the algorithms to be described in the next section, here we explain how 
to optimize a single tensor of the MERA. 

Let if be a Hamiltonian made of nearest neighbor, two-site interactions /i^'^''""^^!, 

if = ^ (5.57) 

r 

For purposes of the optimization below, we choose each term hy^"^^^^ so that it has no 
positive eigenvalues, h^''''^^^^ < 0. [This can be achieved with the simple replacement 
fiV,r+i\ _i. f^[r,T+i\ _ x^^j^ where Amax is the largest eigenvalue of hy''''^^^] 
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Our goal for the time being will be to minimize the energy (Fig. 5.12 i) 



E = tT{HP), (5.58) 

where P is a projector onto the xy-dimensional subspace Yu G by modifying only one 
of the tensors of the MERA. The optimization of a disentangler u is very similar to that 
of an isometry w, and we can focus on describing the latter in more detail. 

Suppose then that, given a MERA, we want to optimize an isometry w while keeping 
the rest of the tensors fixed. The cost function E is quadratic in w (more specifically, it 
depends bi-linearly on w and w'^), 

E{w) = tr(^ wMsW^N,) + ci, (5.59) 

s 

where Mg and Ns are two sets of matrices and ci is a constant (that originates in all the 



Hamiltonian terms of Eq. 5.57 outside the future causal cone of w). Unfortunately there 



is no known algorithm to solve a quadratic problem subject to the additional isometric 



constraint of Eq. 2.2 One can, however, attempt several approximate strategies. Here 



we describe an iterative approach based on linearizing the cost function E{w). 

In this approach, we temporarily regard w and as independent tensors, and optimize 
w while keeping w"^ fixed. The cost function reads, up to the irrelevant constant, simply 

E'-iw) = tr(^^;T^), T^-J^ Ms w^N,, (5.60) 

s 

where we call the matrix the environment of the isometry w and we treat it as if it 
was indepedent of w. E*{w) is then minimized by the choice w = —WV\ where V and W 
are the unitary transformations in the singular value decomposition of the environment, 
= VSW\ 

minE*{w) = mm{wVSW^) = -tr(^) = - V (5.61) 

a 

(here Sq, > are the singular values of T^^,). 

Accordingly, given an initial isometry w, the optimization is performed by iterating the 
following four steps g^ne times: 



LI. Compute the environment T^; with the newest version of (as explained below. 



see also Fig 5.13 ). 
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H 



(ii) Isometry Environment: Y^^ 




Figure 5.12: (i) The energy of a MERA, defined E = tr(ifP), is represented explicitly 
in terms of a tensor network. The removal of an isometry w from this network gives (ii) 
the environment for w (and similarly for disentanglers u). By construction we have 
that E = tr(wT^). 
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L2. Compute the singular value decomposition = VSW'^. 
L3. Compute the new isometry w' = —WV^. 
L4. Replace w"^ with w'"^ . 



Enviromentsof an isometry: 




Figure 5.13: Tensor network corresponding to the 6 different contributions to the 
environment = X]^=i '^w isometry w. Notice that at each iteration of L1-L4 

we need to recompute each TJ„ since it depends on the updated Nevertheless, the 
Hamiltonian term and density matrix that appears in remain the same throughout 
the optimization and only need to be computed once. 

The environment of an isometry w (at layer r) can be decomposed as the sum of 6 
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contributions {i = 1, ■ ■ ■ ,6), each one expressed as a tensor network that involves 
neighboring isometric tensors of the same layer r (disent anglers and isometries) as well as 
one Hamiltonian term h^^f^^^^ and one density matrix pt '^^\ see Fig. 5.13, The two-site 
Hamiltonian term collects contributions from all the Hamiltonian terms in Eq. 



5.57 included in the future causal cone of the sites r, r + 1 of C^-i and is computed with 
the help of the ascending superoperator A. Similarly, the two-site density matrix ' 
is computed with the help of the descending superoperator V. The computation of 
and pr which only needs to be performed once during the optimization of w, has a 

cost O(x^logiV). 

On the other hand, once we have /iJ.T]^^^ and pr computing T^^ has a cost 0{x^) and 
needs to be repeated at each iteration of the steps L1-L4, with a total cost O^x^q^nc)- In 
actual MERA simulations we find that the cost function E^j typically drops very close to 
the eventual minimum already after a small number of iterations q^,^^ of the order of 10. 

The optimization of a disentangler u is achieved analogously, but in this case the environ- 



ment Tu decomposes into three contributions (i = 1, 2, 3), see Fig. 5.14, The required 
Hamiltonian terms and density matrices can be computed at a cost O(x^logA^), while 
the optimization of u following steps L1-L4 has a cost 0{x^qone)- 



5.4 Optimization of the MERA 

In this section we explain a simple algorithm to optimize the MERA so that it minimizes 



the energy of a local Hamiltonian of the form Eq. 5.57 We first describe the algorithm for 
a generic system, and then discuss a number of specialized variations. These are directed 
to exploit translation invariance, scale invariance and to simulate systems where there is 
a finite range of correlations. 



5.4.1 The algorithm 



The basic idea of the algorithm is to attempt to minimize the cost function of Eq. 5.58 by 
sequentially optimizing individual tensors of the MERA, where each tensor is optimized 
as explained in the previous section. 



By choosing to sweep the MERA in an organized way, we are able to update all its 0{N) 
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Environments of a disentangler: 




Figure 5.14: Tensor networks corresponding to the 3 different contributions to the 
environment T„ = "^2^=1 disentangler u. Notice that at each iteration of L1-L4 

we need to recompute each since it depends on the updated . Nevertheless, the 
Hamiltonian term and density matrix that appears in T„ remain the same throughout 
the optimization and only need to be computed once. 
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tensors once with cost 0{x^N). Here we describe a bottom-top approach where the MERA 
is updated layer by layer, starting with the bottom layer r = 1 and progressing upwards 
all the way to the top layer (top-bottom and combined approaches are also possible). 



Given a starting MERA and the Hamiltonian of Eq. 5.57 a bottom-top sweep is organized 
as follows: 



Al. Compute all two- site density matrices pr for all layers r and sites r G 



Starting from the lowest layer and for growing values of r = 1, 2, ■ ■ ■ , T — 1, repeat the 
following two steps: 

A2. Update all disentanglers u and isometrics w of layer r. 
A3. Compute all two-site Hamiltonian terms /ir'*^^^^ for layer r. 



Then, finally. 



A4. Update the top tensor of the MERA. 



In step Al, we compute all nearest neighbor reduced density matrices pr''^'^^^ for all the 
lattices Cr, so that they can be used in step A2. We first compute the density matrix for 



the two sites of Ct-i as explained in Fig. 5.5 Then we use the descending superoperator 



V to compute the 6 possible nearest neighbor, two-site density matrices of Ct-2- More 
generally, given all the relevant density matrices of Cr, we use V to obtain all the relevant 
density matrices of £r-i- In this way, the number of operations is proportional to the 
number of computed density matrices, namely 0{N), and the total cost is 0{x^N). 

Step A2 breaks into a sequence of single-tensor optimizations that sweeps a given layer r 
of the MERA. Each individual optimization in that layer is performed as explained in the 
previous section. Note that in order to optimize, say, an isometry w, we build its envi- 
ronment Tw by using (i) the density matrices computed in step Al; (ii) the Hamiltonian 
terms that were either given at the start for r = 1 or have been computed in A3 for r > 1; 
(iii) the neighboring disentanglers and isometrics within the layer r. We can proceed, for 
instance, by updating all disentanglers of the layer from left to right, and then update all 
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the isometries. This can be repeated a number q^^^ of times until the cost function does 
not change significantly. 

In step A3, the new disentanglers and isometries of layer r are used to build the ascending 
superoperator which we then apply to the Hamiltonian terms of layer r — 1 to compute 



the Hamiltonian terms for layer r. As explained after Eq. |5.16[ each Hamiltonian term 
in layer r is built from three contributions from layer r — 1 . 

In step A4, the optimized top tensor corresponds to the xt eigenvectors with smaller 
energy eigenvalues of the Hamiltonian of the two-site lattice Ct-i, obtained by exact 
diagonalization. 

The overall optimization of the MERA consists of iterating steps A1-A4 until some pre- 
established degree of convergence in the energy E is achieved. Suppose this occurs after 
giter iterations. Then the cost of the optimization scales as 0(x'^A^5'onc5'iay5'itor)- We observe 
that it is often convenient to keep gone and qi^y relatively small (say between 1 and 5), 
since it is not worth spending much effort optimizing a single tensor /layer that will have 
to be optimized again later on with a modified cost function. 

5.4.2 Translation invariant systems 

When the Hamiltonian H is invariant under translations, we can use a translation invari- 
ant MERA 

In this case, each layer r is characterized by a disentangler Ur and an isometry Wr- In 
addition, on each lattice Cr we have one two-site hamiltonian hr and one average density 
matrix p^. Then a bottom-top sweep of the MERA breaks into the steps A1-A4 for the 
inhomogeneous case above, but with the following simplifications: 

In step Al, we compute Pr-i from using the average descending superoperator V of 



Eq. 5.34 



Pr ^ Pr-1 = ViPr). (5.62) 

Then the whole sequence {pT-ir ' ' iPiiPq}i with T ^ logA^, is computed with cost 



translation invariant MERA, characterized by one disentangler and one isometry at each layer, 
need not represent a translation invariant state |\['). Hence the need to consider the average density 
matrix p 
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In step A2, the minimization of the energy E by optimizing, say, the isometry Wr is no 
longer a quadratic problem (since a larger power of Wr appears now in the cost function). 
Nevertheless, we still linearize the cost function E and optimize Wr according to the 
steps L1-L4 of the previous section. Namely, we build the environment T^^, (which now 
contains copies of Wr and wl, all of them treated as frozen), compute its singular value 
decomposition to build the optimal w'^, and then replace Wr and wl with w'^ and w'^ in 
the tensor network for before starting the next iteration. 

In step A3, the new hamiltonian term hr is obtained from hr-i using the average ascending 
superoperator A, 

hr-l ^ hr = A{hr-l). (5.63) 

Step A4 proceeds as in the inhomogeneous case. 

The overall cost of optimizing the MERA is in this case 0(x^ log(A^)gonc'?lay5'itcr)• 
5.4.3 Scale invariant systems 

Given the Hamiltonian H for an infinite lattice at a quantum critical point, where we 
expect the system to be invariant under rescaling, we can use a scale invariant MERA 
to represent its ground state |Vid07l IVidOSl lEVlObl lEVlOal IGMF081 IPEV09j . A scale 
invariant MERA is also relevant in the context of topological order in the infra-red limit 
of the RG flow [AVOSj IKRV09j , both for finite and infinite systems; we will not consider 
such systems here. 

Let us assume that all disentanglers and isometrics are copies of a unique pair {u,w). 
Then, as explained in Ref. |PEV09j . the optimization algorithm can be specialized to 
take advantage of scale invariance as follows: 

In step Al, we apply sparse diagonalization techniques to compute the fixed point density 
matrix p from the superoperator S*. This amount to applying S* a number of times and 
therefore can be accomplished with cost 0{x^). 

In step A2, the environment for e.g. the isometry w, T^, is computed as a weighted sum 
of environments for different layers t = 1,2, ■ ■ ■ . In a translation invariant MERA the 
environment for layer r is a function f{ur,Wr,pT,hr-i) of the pair {ur,'Wr), the density 
matrix pr and the Hamiltonian term hr-i (specifically, / is the sum of the diagrams in 
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Fig. 5.13). A scale invariant MERA corresponds to the replacements 

{Ur,Wr) {U,W), Pr P, (5.64) 

SO that only hr-i retains dependence on r. We then choose the average environment 

oo ^ 

'^^ = ^^f('^^^^P^^r-i), (5.65) 



3- 

T=l 



where the weight 1/3"^ reflects the fact that for each isometry at layer r there are 3 
isometrics at layer r — 1. Using linearity of the / in its fourth argument we arrive at 

oo 

T^ = f{u,w,p,h), /i = y'— (5.66) 

r=l 

In practice, only a few terms of the expansion of h (say r = 1, 2, 3, 4) seem to be necessary. 
Given T^„, the optimization proceeds as usual with a singular value decomposition. 

Steps A3 and A4 are not necessary. 

That is, the algorithm to minimize the expected value of H consists simply in iterating 
the following two steps: 



Sclnvl. Given a pair {u,w), compute a pair (p, h). 

Sclnv2. Given the pairs {u,w) and (p, h), update the pair {u,w). 



In practical simulations it is convenient to include a few (say one or two) transitional layers 
at the bottom of the MERA, each one characterized by a different pair (m^, Wr)- In this 
way the bond dimension x of the MERA can be made independent of the dimension d of 
the sites of C. These transitional layers are optimized using the algorithm for translation 
invariant systems. 

5.4.4 Finite range of correlations 

A third variation of the basic algorithm consists in setting the number T' of layers in the 
MERA to a value smaller than its usual one T ^ logg N (in such a way that the number 
Nt' of sites on the top lattice Ct' may still be quite large) and to consider a state of 
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Finite range MERA: 




Figure 5.15: A finite correlation range MERA with T' = 2 layers is used to represent 
a state of = 36 sites. Since it lacks the uppermost layers, only sites within a finite 
distance or range C ~ 3^ of each other may be correlated. More precisely, only pairs of 
sites (ri, whose past casual cones intersect can be correlated, as it is the case with the 
pair of sites (6, 14) but not with (14, 26), for which we have (ol^^lot^^l) = (o[^^])(o[26]). 



the lattice C such that after T' coarse-graining transformations it has become a product 
state, 

l^o) ^ l^i) ^ > I^T') , (5.67) 



where I^^t') = lO)*^^^'. For instance. Fig. 5.15 shows a MERA for = 36 and T' 



The four top tensors of this MERA are of type (0,3), where the lack of upper index 
indicates that the top lattice Ct' is in a product state of its Nt' = 4 sites. 

We refer to this ansatz as the finite range MERA, since it is such that correlations in |^) 
are restricted to a finite range (, roughly ( ^ 3^ sites, in the sense that regions separated 
by more than ( sites display no correlations. This is due to the fact that the past causal 



cones of distance regions of C have zero intersection, see Fig. 5.15 



Given a ground state |\E') with a finite correlation length ^, the finite range MERA with 
C ~ ^ turns out to be a better option to represent than the standard MERA with T ^ 
log3 layers, in that it offers a more compact description and the cost of the simulations 
is also lower since there are less tensors to be optimized. The algorithm is adapted in a 
straightforward way. The only significant difference is that the top isometrics, being of 
type (0,3), do not require any density matrix in their optimization (their environment is 
only a function of neighboring disentanglers and isometries, and of Hamiltonian terms). 

A clear advantage of the finite range MERA is in a translation invariant system, where 
the cost of a simulation with range ( = 3^' is 0(x^ logs C); that is, independent of A^. This 
allows us to take the limit of an infinite system. We find that, given a translation invariant 
Hamiltonian H = Yl^=i h^^'^~^^\ where is the same for all r E C, the optimization 
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of a finite range MERA will lead to the same collection of optimal disentanglers and 
isometrics {{ui,wi), {u2, W2), ■ ■ ■ , (mt', wt')}, for different lattice sizes A^, A^', A^" ■ ■ ■ larger 
than (. This is due to the existence of disconnected causal cones, which imply that the 
cost functions for the optimization are not sensitive to the total system size provided it 
is larger than (. As a result, {(ui, Wi), {u2, W2), ■ ■ ■ , {ut', wt')} can be used to define not 
just one but a whole collection of states |^'(A^)), |^(A^'))> for lattices of 

different sizes A^, A^', A^", ■ ■ ■ , such that they all have the same two-site density matrix p 
and therefore also the same expected value of the energy per link. 



In particular, we can use the finite range MERA algorithm to obtain an upper bond for 
the ground state energy of an infinite system, even though only T' pairs {ut-,Wt-) are 
optimized. 



5.5 Benchmark calculations for ID systems 

In order to benchmark the algorithms of the previous section, we have analyzed zero 
temperature, low energy properties of a number of ID quantum spin systems. Specifically, 
we have considered the Ising model |Pfe70i IBG85] , the 3-state Potts model |SP81] , the 
XX model |LSM61] and the Heisenberg models |Bax82] , with Hamiltonians 



(^(A^)| h \^{N)) = {^{N')\ h \^{N')) 



(5.68) 




(5.69) 



r 




(5.70) 




(5.71) 



r 




■r 



(5.72) 
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Figure 5.16: Top: The energy error of the MERA approximations to the ground-state 
of the infinite Ising model, as compared against exact analytic values, is plotted both for 
different transverse magnetic field strengths and different values of the MERA refinement 
parameter x- The finite correlation range algorithm (with at most T = 5 levels) was used 
for non-critical ground states, whilst the scale invariant MERA was used for simulations 
at the critical point. It is seen that representing the ground-state is most computationally 
demanding at the critical point, although even at criticality the MERA approximates the 
ground-state to between 5 digits of accuracy (x = 4) and 10 digits of accuracy (x = 22). 
Bottom: Scale-invariant MERA are used to compute the ground-states of infinite, critical. 



ID spin chains of Eqs. 5.69 5.72 for several values of x- instances one observes a 



roughly exponential convergence in energy over a wide range of values for x as indicated 
by trend lines (dashed). Energy errors for Ising, XX and Heisenberg models are taken 
relative to the analytic values for ground energy whilst energy errors presented for the 
Potts model are taken relative to the energy of a x = 22 simulation. 
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where a^-, Uy and are the spin 1/2 Pauli matrices and M^^i-i Mx^2 and are the 
matrices 

/ 2 \ 

M, = 0-10 , (5.73) 
Vo -1/ 

/ 1 \ 



xA 



/ 1 \ 

1 
[o 1 0/ 



(5.74) 



1 

VI 0/ 

We assume periodic boundary conditions in all instances and use a translation invariant 
MERA to represent an approximation to the ground state and, in some models, also 
the first excited state. For Ising and Potts models the parameter A is the strength of 
the transverse magnetic field applied along the z-axis, with Ac = 1 corresponding to 
a quantum phase transition. Both the XX model and Heisenberg model are quantum 
critical as written. 



Fig. 5.16 shows the accuracy obtained for ground-state energies of the above models in 
the limit of an infinite chain, as a function of the refinement parameter x- Simulations 
were performed with either the finite correlation range algorithm (for the non-critical 
Ising) or the scale invariant algorithm (for critical systems). In all cases one observes 
roughly exponential convergence to the exact energy with increasing x- For any fixed 
value of X; the MERA consistently yields more accuracy for some models than for others. 
For the Ising model, the cheapest simulation considered {x = 4) produced 5 digits of 
accuracy, whilst the most computationally expensive simulation {x = 22) produced 10 
digits of accuracy. The time taken for the MERA to converge, running on a 3GHz dual- 
core desktop PC with 8Gb of RAM, is approximately a few minutes/hours/days/weeks 
for X = 4, 8, 16, 22 respectively. We stress that these simulations were performed on single 
desktop computers; a parallel implementation of the code running on a computer cluster 
might bring significantly larger values of x within computational reach. 



Fig. 5.17 demonstrates the ability of the MERA to reproduce finite size effects. It 



shows the transverse magnetization as a function of the transverse magnetic field for 
several system sizes. The results smoothly interpolate between those for small system 
sizes and those for an infinite chain, and match the available exact solutions. On the 
other hand, the MERA can also be used to explore the phase diagram of a system. 



Fig. 5.18 shows the spontaneous magnetization, which is the system's order parameter. 



for a ID chain of = 162 sites for both Ising and Potts models, where has been chosen 
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Figure 5.17: The transverse magnetization (ctz) for Ising and ^ (M^) for Potts, is plotted 
for translation invariant chains of several sizes A^. Top: For the Ising model, the magneti- 
zation given from x = 8 MERA matches those from exact diagonalization for small system 
sizes (A^ = 6, 18), whilst the magnetisation from the = 54 MERA approximates that 
from the thermodynamic limit (known analytically). Bottom: Equivalent magnetisations 
for the Potts model, here computed with a x = 12 MERA. Simulations with larger N 
systems show little change from the = 54 data, again indicating that A^ = 54 is already 
close to the thermodynamic limit. 
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Figure 5.18: Top: Spontaneous magnetization (cr^.) computed with a x = 8 MERA for 
a periodic Ising system of = 162 sites. The results closely approximate the analytic 
values of magnetization known for the thermodynamic limit. A fit of the data near the 
critical point yeilds a critical exponent /3mera = 0.1243, with the exact exponent known 
as /3ex = 1/8. Bottom: An equivalent phase portrait of the Potts model, here with 
spontaneous magnetization ^ (M^. i + M^, 2), is computed with a x = 14 MERA and is 
plotted with a fit of the data near the critical point. The fit yields a critical exponent 
/Smera = 0.105 with the exact exponent known to be /3ex = 1/9. 
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large enough that the results under consideration do not change singnificantly with the 
system size (thermodynamic limit). A fit for the critical exponent of the Ising model 
gives /Smera = 0.1243 whilst the fit for the Potts model produces /3mera = 0.105. These 
values are within less than 1% and 6% of the exact exponents /3 = 1/8 and /3 = 1/9 
for the Ising and Potts models respectively. Obtaining an accurate value for this critical 
exponent through a fit of the data near the critical point is difficult due to the steepness 
of the curve near the critical point. Through an alternative method involving the scaling 



super-operator S (Sect. 5.2), more accurate critical exponents can be obtained, as shall 



be demonstarted in Chapter [6] 

The previous results refer to local observables. Let us now consider correlators. A scale- 
invariant MERA, useful for the representation of critical systems, gives polynomial cor- 



relators at all length scales, as shown in Fig. 5.19 for the critical Ising and Potts models. 



Note that Fig. 5.19 displays the correlators that are most convenient to compute (as per 



Fig. |5.10[ ). These occur at distances ci = 3'^ for g = 1, 2, 3 . . . and are evaluated with 



cost 0{x^). Evaluation of arbitrary correlators is possible (see Fig. 5.8) but its cost is 
several orders of x more expensive. The precision with which correlators are obtained is 
remarkable. A x = 22 MERA for the Ising model gives {ctPctx'^'^) correlators, at distances 
up to (i = 10^ sites, accurate to within 0.6% of exact correlators. Critical exponents t] 
are obtained through a fit of the form C{r,r + d) oc with C as correlators of x,y or 
z magnetization. For the Ising model, the exponents for x, y and z magnetizations are 
obtained with less than 0.02% error each. For the Potts model exponents are obtained 
with less than 0.04% and 0.08% error for x and z magnetization respectively. 

Finally, we also demonstrate the ability of MERA to investigate low-energy excited states 



by computing the energy gap in the Ising and Potts models. Fig. 5.20 shows that the gap 
grows linearly with the magnetic field A and independent of in the disordered phase 
A >> Ac, whilst at criticality it closes as 1/A^. Even a relatively cheap x = 8 MERA 
reproduces the known critical energy gaps to within 0.2% for systems as large as = 162 
sites. The expected value of arbitrary local observables (besides the energy) can also be 
easily evaluated for the excited state. 
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Figure 5.19: Two-point correlators for infinite ID Ising and Potts chains at criticality 
(A = 1), as computed with x = 22 scale- invariant MERA. Correlators for the Ising model 
are compared against analytic solutions |Pfe70l IBG85] whilst those for the Potts model are 
plotted against the polynomial decay predicted from CFT |PDF S97]. A scale-invariant 
MERA produces polynomial decay of correlators at all length scales; a fit of the form 
{cTx^ (Jx'^'^ ) oc d~^'^ for Ising correlators generated by the MERA gives the decay exponent 
rj^ = 0.24996, close to the known analytic value 1/4 and similarly for the fits on other 
correlators. Indeed the MERA here reproduces exact {cTx^cTx'^^'^) correlators for the Ising 
model at a distance up to (i = 10^ sites within 0.6% accuracy. Critcal exponents for the 
Potts are also reproduced very accurately. 
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Figure 5.20: Top: A % = 8 MERA is used to compute the energy gap AE (the energy 
difference between the ground and 1*^* excited state) of the Ising chains as a function of 
the transverse magnetic field. The gap computed with MERA for = 6, 18 sites is in 
good agreement with that computed through exact diagonahzation of the system. Inset: 
Crosses show analytic values of energy-gaps at the critical point for N = {6, 18, 54, 162}. 
Even for the largest system considered, = 162, the gap computed with MERA 
^-Emera = 9.67 X 10~^ compares well with the exact value AE^^ = 9.69 x 10^'^. Bot- 
tom: Equivalent data for the Potts Model where simulations have been performed with a 
X = 14 MERA to account for the increased computational difficulty of this model. 
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5.6 Conclusions 



In this Chapter we have provided a rather self-contained description of an algorithm to 



explore low energy properties of lattice models (Sect. 5.3 5.4), and benchmark calculations 



addressing ID quantum spin chains (Sect. 5.5). 



Many of the features of the MERA algorithm highlighted by the present results were 
also observed in systems of free fermions and free bosons in Chapters |3] and |4| These 
include (i) the ability to consider arbitrarily large systems, (ii) the ability to compute the 
low-energy subspace of a Hamiltonian, (iii) the ability to disentangle non-critical systems 
completely, (iv) the ability to find a scale-invariant representation of critical systems 
and finally (v) the reproduction of accurate polynomial correlators for critical systems. 
However, the algorithms of Chapters [3] and |4] exploit the formalism of Gaussian states 
that is characteristic of free fermions and bosons and cannot be easily generalized to 
interacting systems. Instead, the algorithms discussed in this Chapter can be used to 
address arbitrary lattice models with local Hamiltonians. 

An alternative method to optimise the MERA is with a time evolution algorithm as 
described in Ref. |RMV08j . The time evolution algorithm has a clear advantage: it 
can be used both to compute the ground state of a local Hamiltonian (by simulating an 
evolution in imaginary time) and to study lattice dynamics (by simulating an evolution 
in real time). We find, however, that the algorithm described in the present Chapter is 
a better choice when it comes to computing ground states. On the one hand, the time- 
evolution algorithm has a time step 6t that needs to be sequentially reduced in order to 
diminish the error in the Suzuki- Trotter decomposition of the (euclidean) time evolution 
operator. In the present algorithm, convergence is faster and there is no need to fine tune 
a time step St. In addition, the present algorithm allows to compute not only the ground 
state but also low energy excited states. It is unclear how to use the time evolution 
algorithm to achieve the same. 

The benchmark calculations presented in this manuscript refer to ID systems. For such 
systems, however, DMRG |Whi92t IWhi93] already offers an extraordinarily successful 
approach. The strength of entanglement renormalization and the MERA relies on the fact 
that the present algorithms can also address large 2D lattices, as is discussed Chapter [7} 
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Chapter 6 

Entanglement renormalization and quantum 
criticality 



6.1 Introduction 

The study of quantum critical phenomena through real-space renormalization group (RG) 
techniques |Wil75[ IWhi92[ IWhi93j has traditionally been obstructed by the accumulation, 



over successive RG transformations, of short-range entanglement across block boundaries. 
Entanglement renormalization was proposed as a technique to address this problem. By 
removing short-range entanglement at each iteration of the RG transformation, not only 
can arbitrarily large lattice systems be considered, but the scale invariance characteristic 
of critical phenomena is also seen to be restored. 

In this Chapter we explain how to use the MERA to investigate scale invariant systems. 
It has recently been shown that the scale invariant MERA can represent the infra-red 
limit of topologically ordered phases |AV08l IKRV09] , here we focus instead on its use at 
quantum criticality. We present the following results: (i) given a critical Hamiltonian, an 
adaptation of the algorithm of the previous Chapter to compute a scale invariant MERA 
for its ground state; then, starting from a scale invariant MERA, (ii) a procedure to 
identify the scaling operators/dimensions of the theory and (iii) a closed expression for 
two-point and three-point correlators; (iv) a connection between the MERA and conformal 
field theory, which can be used to readily identify the continuum limit of a critical lattice 
model; finally (v) benchmark calculations for the Ising and Potts models. 

We note that result (ii) was already discussed by Giovannetti, Montangero and Fazio in 
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Ref. |GMF08j using the binary MERA of Ref. |Vid08] . Our derivations, as with those of 
the previous Chapter, are conducted using the ternary MERA, in terms of which results 
(iii)-(iv) acquire a simple form. 

We start by considering a finite ID lattice C made of sites, each one described by a 
vector space V of dimension x- The (ternary) MERA is a tensor network that serves 



as an ansatz for pure states |\E') G V®^ of the lattice, see Fig. 6.1 Its tensors, the 
disentanglers and isometries, are organized in T log3 A^ layers, each one implementing 
a RG transformation. Such transformations produce a sequence of lattices, 

Co ^ Ci ^ ■■■ ^ Ct, Co = C, (6.1) 

where lattice Cr+i is a coarse-graining of lattice Cr, and the top lattice Ct is sufficiently 
small to allow exact numerical computations. Let o denote a local observable supported on 
two contiguous sites of £, and let pr be the density matrix that describes the state of the 
system on two contiguous sites of Ct- Then the ascending and descending superoperators 
Ar and 

Or = Ar{Or-l), Pr-1 = T^t{Pt), (6.2) 

generate a sequence of operators and density matrices 

oo ^ oi ^ • • ■ ^ ot, oo = o, (6.3) 

Po <- Pi ^ ■ ■ ■ <- PT, Po = P, (6.4) 



where Or and p^ are supported on two contiguous sites of the lattice Cr- Eq. (6.3) allows 



us to monitor how the local operator o transforms under successive RG transformations. 



whereas its expected value (o) = tr(po) can be evaluated by computing p in Eq. (6.4) 



6.2 RG fixed point 

The scale invariant MERA corresponds to the limit of infinitely many layers, T — )■ oo, 
and to choosing the disentanglers and isometries in all layers to be copies of a unique pair 
u and w [VidOTj IVid08] . In this case we refer to the ascending superoperator Ar, which 



no longer depends on r, as the scaling superoperator S (see Fig. 6.1), and to its dual 
as S*. Notice that iS is a fixed-point RG map. Then, as customary in RG analysis 
|Car96[ IPDFS97] , the scaling operators (pa and scaling dimensions of the theory. 



5(0q) = \a(pa, Aa = - logg K, 



(6.5) 



6.3 Correlators 
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are obtained by diagonalizing this map, 

S{») = ^ Aa0„tr(0a»), tr(0a0/3) = (6.6) 



6.6 



was first discussed 



where (pa are the eigenvectors of the dual S*, S*{(f)a) = ^a<Pa- Eq- 
in Ref. |GMF08j by Giovannetti, Montangero and Fazio. It formalizes a previous obser- 
vation (see Eq. 5 of Ref. |Vid08j ) that the scale invariant MERA displays polynomial 
correlations. By construction, S is unital, S{I) = I, so that the identity operator I in V®^ 
is a scahng operator with eigenvalue Aj = 1; and contractive, meaning |Aq,| < 1 |BR79] . 
Here we will assume, as it is the case in the examples below, that only the identity opera- 
tor I has eigenvalue A = 1. Then the operator p = t is a density matrix that corresponds 
to the unique fixed point of S*, S*{p) = p, and since 

hm ( S*o- - -oS*, ){pt)=P (6.7) 

T— >-oo V ' 

T times 

for any starting px, it follows that p is the state of any pair of contiguous sites of C 
(consistent with scale invariance, p is also the state of any pair of contiguous sites of Cr 
for any finite r). The computation of the expected value of the local observable o is then 
straightforward, 

(o)=tr(/3o), (6.8) 
which for the scaling operators reduces to {(pa) = ^ai- 



6.3 Correlators 



Let us now diagonalize the one-site scaling superoperator S^^^ of Fig. 6^ 



(6.9) 



where the scaling dimensions Aa ■* = — log3 Aa ■* coincide with The correlator for 



two scaling operators (pa^ and (p^^^' placed on contiguous sites reads 

Ca, ^ U'\l)4\0)) = tr((0« ® 0(^))p). 



Suppose now that (pa^ and (f))^' are placed in two special sites x, y as in Fig. 
X — y is such that \rxy\ = 3'' for q = 1,2, ■ ■ ■ . Then after q = logg \r, 



xy 



xy 



6.2 



(6.10) 



where 



iterations 



^Our numerics show that the lowest ua scaling dimensions fulfill ^ ~ ~ A^^"^, where ha grows 
with X- 
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Figure 6.1: (i) Two lowest rows of disentanglers u and isometries w of the ternary MERA. 
They map the original infinite lattice Co = C into increasingly coarse-grained lattices Ci 
and £2- Notice that three sites of Cr-i become one site of Cr, hence the use of base 
3 logarithm throughout the Chapter. {ii)-{w) Under the coarse-graining transformation 
defined by the MERA, two-site operators supported on three different pairs of sites of 
Cr-i become supported on the same pair of sites of Cr- (v) Accordingly, the scaling 
superoperator S is the average of three contributions, each of which (and thus also their 
average) is unital and contractive thanks to the isometric character of u and w |Vid08] . 



6.4 Conformal field theory 
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of tlie RG transformation, 0q ^ and 0^^'' become first neiglibors again. Notice tliat eacli 
iteration contributes a factor Aa^A^^''. Using tlie identity a^°^^ = 6^°^" g^d 

(XW\W)^oSi3\r.y\ - \r |log3(AL'>A(,i') _ I \-Ai'^-Af^ 
K^a J — \'xy\ — I'xyl 

and obtain a closed expression for two-point correlators, 

For three-point correlators we define the constants 

^ AW + A(^)-A« (6.12) 
C^p, ^ 2^-tr((0W®4^)®</)(i))p(3)) (6.13) 

where the trace corresponds to the correlator on three consecutive sites and p^^^ is obtained 
from p. For \rxy \ = \ryz\ = \rxz\/2 = 3'^, analogous manipulations lead to 




(6.14) 



Figure 6.2: (i) One-site operators on special sites are coarse-grained into one-site opera- 
tors, (a) Scaling superoperator for one-site operators, {iii) In computing correlators on 
specific sites x and y (or x, y and z), one-site operators are coarse-grained individually 
according to S^^^ until they become nearest neighbors (which in this case occurs at lattice 
C2,q = 2). 



6.4 Conformal field theory 



The continuous limit of a quantum criticial lattice system (scale invariant case) corre- 
sponds to a conformal field theory (CFT) |Car96t IPDFS97] . A CFT contains an infinite 
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set of quasi-primary fields 0^^^, with scaling dimensions A^^"^. The correlators involving 



two or three quasi-primary fields have expressions analogous to Eqs. |6.11| and |6.14[ and 
the (symmetric) coefficients C^^'^ for three-point correlators coincide with those in the 
so-called operator product expansion (OPE). Moreover, quasi-primary fields are organized 
in conformal towers corresponding to irreducible representations of the Virasoro algebra. 
Each tower contains one primary field 0^ at the top, with conformal dimensions {t, t) 
(such that its scaling dimension is = t + t), and its infinitely many descendants, which 
are quasi-primary fields with scaling dimension A = A^ + n for some integer n > 1. 

A CFT is completely specified by its symmetries once the following conformal data has 
been provided: (i) the central charge c, (ii) a complete list of primary fields with their 
conformal dimensions and (iii) the OPE for these primary fields. For instance, the Ising 
CFT in 1+1 dimensions has central charge c = 1/2, three primary fields identity I, spin 
a and energy e with conformal dimensions (0,0), {j^, j^) and (|, ^), and OPE coefficients 

/^CFT r ^CFT ^ ^CFT ^CFT ^CFT ri ( (1 ^ K\ 



The present analysis readily suggests a correspondence between the scaling operators 0^ 
of the scale invariant MERA, defined on a lattice, and the quasi-primary fields 0^^^ of 
a CFT, defined in the continuum. Together with the algorithm described below, this 
correspondence grants us numerical access, given a critical Hamiltonian H on the lattice, 
to most of the conformal data of the underlying CFT, namely to scaling dimensions and 
OPE coefficients. The central charge c can also be obtained from the von Neumann 
entropy S{p) = — tr(plog2p), which for a block of L sites scales, up to some additive 
constant, as 5 = f logsL jVLRKOSl ICCOi] . We then have S{p) - S{p^^^) = f (log2 2 - 
log2 1) = f , or simply 

c = 3(5(p)-5(p«)). (6.16) 



6.5 Algorithm for scale invariant MERA 

Given a critical Hamiltonian H for an infinite lattice, we obtain a scale invariant MERA 
for its ground state by adapting the general strategy discussed in Chapter |5} Recall 
that tensors (disentanglers u and isometrics w) are optimized so as to minimize the energy 
E = After linearization this reads 

E = ti{uT^) + ki = tr{wT^) + k2, (6.17) 



6.6 Benchmark results 
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where T„ and T^„ are known as environments and fci, ^2 are two irrelevant constants. In 
the translation invariant case the environment for, say, an isometry w at layer r of the 
MERA, = f{ur,Wr,pT,hr-i), is a function of the disentangler Ur and isometry Wr 
of that layer, a two-site density matrix p^- and a two-site Hamiltonian term hr-i- In the 
present case, we replace the above with the unique pair {u,w), the fixed-point density 
matrix p, and an average Hamiltonian h = ^^h^/y ^ where the weights 1/3"^ account 
for the relative number of tensors in different layers of the MERA. Then, starting from 
some initial pair (m, w) and the critical Hamiltonian H made of two-body terms h, the 
following steps are repeated until convergence: 

1. Given the latest {u,w), compute {p,h). 

2. Given {u,w,p,h), update the pair {u,w). 



In step Al, the scaling superoperator S is built as indicated in Fig. |6.1[ We compute the 
fixed-point density matrix p by sparse diagonalization of S, and the average Hamiltonian 
h by using hr = S{hr-i), ho = Step A2 is decomposed into a sequence of alternating 
optimizations for u and w as in the generic algorithm of Chapter |5} where each tensor is 
updated by computing a singular value decomposition of its environment. 



6.6 Benchmark results 

We illustrate the above ideas and the performance of the algorithm by considering the 
Ising and 3-level Potts quantum critical models in ID, 

HMn, = (6-18) 

r 
r 

where cr^ and ax are Pauli matrices, and 
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^In practice we only compute the first k terms (fc « 2, 3) of tire expansion h = Iiq + hi/3 + h2/9 + ■ ■ ■ ■ 
Tliis average is only needed when H contains operators that are irrelevant in the RG sense. 
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Figure 6.3: Scaling dimensions Aq, obtained from the spectrum of the scahng superop- 
erator S. Circles indicate primary fields. Left: For the Ising model we can identify the 
scaling dimensions of the three primary fields, the so-called identity I, spin a and energy 
e, together with several of their descendants. Right: The spectrum of S for the 3-level 
Potts model shows some of its primary fields, including its primary fields with multiplicity 
two, namely the spins ai and a2 and the pair Zi and Z2 |PDFS97] . 



Mx^2 = {Mx^iy. Notice that sites have a vector space of dimension d = 2 or d = 3. In 
order to use a scale invariant MERA with x > d, we allow the disentanglers and isometrics 
of the first few (typically one to five) layers to be different from u and w. We iterate steps 
A1-A2 about 1000 times. With a cost per iteration that scales as ^ind using a 3 GHz 
dual core desktop with 8 Gb of RAM, simulations for x = 4, 8, 16, 22 take of the order of 
minutes, hours, days and weeks respectively. The following results correspond to x = 22. 



From Eq. 6.16 we obtain an estimate for the central charge, namely cising = .5007 and 



cpotts = -806, to be compared with the exact results 0.5 and 0.8. Fig. |6.3| shows the 
smallest scaling dimensions of the scaling superoperator S We obtain remarkable 
agreement with those expected from GET, as shown in Table |6.1[ Recall that all the 
critical exponents of the model can be obtained from the scaling dimensions of primary 
fields. For instance, for the Ising model the exponents z/ and r] are u = 2Ao- and r] = , 
whereas the scaling laws express the critical exponents a,/?, 7, 5 in terms of u and 77 
[PDFS97] . Further, the OPE coefficients for primary fields of, say, the critical Ising 



model are computed as follows. The matrix Cas in Eq. 6.10 is diagonal for the scaling 



operators corresponding to I, a and e, which we normalize so that Cap = 5ap- With this 



■'Our numerics show that the lowest ua scaling dimensions fulfill ^ « « A^^"^, where grows 
with X- 



6.7 Discussion 
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Table 6.1: Comparison of scaling dimensions of primary fields of the Ising and Potts 
models calculated using MERA (A(mera x = 22)) with exact results known from CFT 



normalization, we then compute the coefficients Cap-y using Eq. 6.13 We reproduce all 
the values of Eq. 6.15 with errors bounded by 3 x 10"''. 



6.7 Discussion 

In this Chapter we have explained how to compute the ground state of a critical Hamil- 
tonian using the scale invariant MERA and how to extract from it the properties that 
characterize the system at a quantum critical point. Our results, which build upon those 
of Ref. |Vid07t lEVlObi lEVlOai IVid08[ EV08| I(;MF08[ IEV09aj . also unveil a concise con- 
nection between the scale invariant MERA and CFT. This correspondence adds signifi- 
cantly to the conceptual foundations of entanglement renormalization. The scale invariant 
MERA can be regarded as approximately realizing an infinite dimensional representation 
of the Virasoro algebra |Car96t [PDFS97j . The finite value of x effectively implies that only 
a finite number of the quasi-primary fields of the theory can be included in the descrip- 
tion. Fields with small scaling dimension, such as primary fields, are retained foremost. 
As a result, given a Hamiltonian on an infinite lattice, we can numerically evaluate the 
scaling dimensions and OPE of the primary fields of the CFT that describes the contin- 
uum limit of the model. This approach differs in a fundamental way from, and offer an 
alternative to, the long-established techniques of Refs. jCar84l ICar86j . based instead on 
finite size scaling. We conclude by noting that most of our considerations rely on scale 
invariance alone and can be applied to study also critical ground states in 2D systems, as 
is demonstrated next in Chapter [7} 
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Chapter 7 



Entanglement renormalization in two spatial 
dimensions 



7.1 Introduction 

Entanglement renormalization has been proposed as a real-space renormalization group 
(RG) method |Wil75] to study extended quantum systems on a lattice. A highlight of the 
approach is the removal, before the coarse-graining step, of short-range entanglement by 
means of unitary transformations called disentanglers. This prevents the accumulation 
of short-range entanglement over successive RG transformations. Such accumulation is 
the reason why the density matrix renormalization group (DMRG) |Whi92t IWhi93] - an 
extremely powerful technique for lattices in one spatial dimension - breaks down in two 
dimensions, where it can only address small systems. 

The use of disentanglers leads to a real-space RG transformation that can in principle be 
iterated indefinitely, enabling the study of very large systems in a quasi-exact way. This 
RG transformation also leads to the so-called multi-scale entanglement renormalization 
ansatz (MERA) |Vid08j to describe the ground state of the system - or, more generally, 
a low energy sector of its Hilbert space. In a translation invariant lattice made of N sites, 
the cost of simulations grows only as logA^ |EV09aj . In the presence of scale invariance, 
this additional symmetry is naturally incorporated into the MERA and a very concise 
description, independent of the size of the lattice, is obtained in the infrared limit of a 
topological phase |AV08l IKRV09] or at a quantum critical point |Vid07l IVid08| lEVlObl 
lEVlOal [GMF081 IPEV091 lMRGF09j . 
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The basic principles of entanglement renormalization are the same in any number of spatial 
dimensions, however, numerical work with 2D lattices incurs a much larger computational 
cost. In Chapters [3] and |4] we explored the use of ER in 2D lattices of free fermions and 
free bosons. Entanglement Renormalization has also been tested for an Ising model in a 
square lattice of small linear size L < 8 |CDR08j . It must be emphasized, however, that 
the approach of Chapters [3] and |4] relies on the gaussian character of free particles and 
can not be generalised to the interacting case, whereas the results of Ref. |CDR08j were 
obtained by exploiting a significant reduction in computational cost that occurs only for 
small 2D lattices. 



In this Chapter we present an implementation of the MERA that allows us to consider, 
with modest computational resources, 2D systems of arbitrary size, including infinite sys- 
tems. In this way we demonstrate the scalability of entanglement renormalization in two 
spatial dimensions and decisively contribute to establishing the MERA as a competitive 
approach to systematically address 2D lattice models. The key of the present scheme is 
a carefully planned organization of the tensors in the MERA, leading to simulation costs 
that grow as 0(x^^), where x is the dimension of the vector space of an effective site. This 
is drastically smaller than the cost 0(x^*) of the 4-to-l MERA scheme (see Fig. 2.7) that 



was used in Chapters |3] and |4] and in Ref. |CDR08] . We also demonstrate the performance 
of the scheme by analysing the 2D quantum Ising model, for which we obtain accurate 
estimates of the ground state energy and magnetizations, as well as two-point correlators 
(shown to scale polynomially at criticality), the energy gap, and the critical magnetic field 
and beta exponent. Finally, we discuss how the use of disentanglers affects the simulation 
costs, by comparing the MERA with a tree tensor network (TTN) |TEV09j . 



7.2 The MERA on a 2D lattice 

Let us consider a square lattice £0 made oi N = L x L sites, each one described by 
a Hilbert space V of finite dimension d. The proposed 2D MERA is characterized by 



the coarse-graining transformation of Fig. (7.1), where blocks of 3 x 3 sites of lattice 



Co are mapped onto single sites of a coarser lattice £1. This is achieved in three steps: 
first disentanglers u are applied on the four sites located at the corners of four adjacent 
blocks; then disentanglers v are applied at the boundary between two adjacent blocks, 
transforming four sites into two; finally, isometrics w are used to map a block into a single 
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Figure 7.1: Entanglement renormalization scheme for a square lattice. A block of 3 x 3 
sites of lattice £r-i (i) is mapped onto one site of Cr (v). The RG transformation 
involves (ii) applying disentanglers u between the corners of adjacent blocks followed by 
(iii) disentanglers v which act across the sides of adjacent blocks and (iv) isometrics w 
which act within a block. Tensors m, v and w have a varying number of incoming and 



outgoing indices (vi) according to Eq. 7.1 
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effective site. In this way, tensors m, v and w 

yt.Y®4^y®4^ ^t.Y®4^y®2^ y^t . _^ V, (7.1) 

transform the state |\l/o) ^ V®^ of the lattice £o in which we are interested (typically the 
ground state of a local Hamiltonian Hq) into a state |\E'i) G V®^/^ of the effective lattice 
£i through the sequence 

(7.2) 

To understand the role of these tensors, it is useful to think of the state |\I'o) as possessing 
three different kinds of entanglement: short-range entanglement residing at the corners 
of four adjacent blocks, short-range entanglement residing near the boundary shared by 
two blocks, and long-range entanglement. Then the disentanglers u and v are used to 
reduce the amount of short-range entanglement residing near the corners and boundaries 
of the blocks. In other words, in states |\E'q) and |\E'o) increasing amounts of short-range 
entanglement from |\E'o) have been removed. This fact facilitates significantly the job of 
the isometry w, namely to compress into an effective site of £i those degrees of freedom in 
a block that still remain entangled (now mostly through long-range entanglement) with 
degrees of freedom outside the block. Thus, the resulting state |\E'i) still contains the 
long-range entanglement of l^&o), but most of its short-range entanglement is gone. We 
complete the above construction by noticing that a dimensional space V is often too 
small to accommodate all the relevant degrees of freedom left on a block. Accordingly, we 
shall describe the effective sites of Ci with a space of larger dimension %. This dimension 
X determines both the accuracy and cost of the simulations. 



The transformation of Fig. 7.1 can now be applied to lattice £i, producing a coarser 
lattice £2- More generally, if £0 is finite, O(logA^) iterations will produce a sequence of 
lattices {Cq, Ci, C2, ■ ■ ■ , Aop} where the top lattice Aop contains only a small number of 
sites and can be addressed with exact numerical techniques. Thus, given a Hamiltonian 
on £05 we can use the above RG transformation to obtain a sequence of Hamiltonians 
{Hq, Hi, H2, ■ ■■ ,iftop}, then diagonalize H^^^ to find its ground state I^E'top), and finally 
recover the ground state |\E'o) of Hq by reversing all the RG transformations: 

l^top) ^ y 1^2) ^ 1^1) ^ 1^0) • (7.3) 



^The tensors of the MERA are called disentanglers u or isometries w depending on whether they 
are in charge of eliminating short-range entanglement or of mapping a block of sites into a single site 
[Vid07| IVid08|, IEV09a| . This distinction is somewhat arbitrary: one can consider tensors that fulfill the 
two roles simultaneously, such as tensor v in Fig. |7.1[ that we still call disentangler. All these tensors 



must be isometric, that is u^u = I, w^w = J, w^w = I. The hermitian conjugation (f) in Eq. (7.11 
appears for consistency with previous references 
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This is precisely how the MERA is defined. Specifically, the MERA for |\E'o) is a tensor 
network containing (i) a top tensor, that describes |^I/)j^p, and (ii) O(logA^) layers of 
tensors (disentanglers and isometrics), where each layer is used to invert one step of the 



coarse-graining transformation of Fig. 7.1 according to the sequence (7.3). 
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Figure 7.2: Spontaneous and transverse magnetizations (a^) and (a^) as a function of 
the applied magnetic field A and for different lattice sizes L. Results for small systems 
correspond to exact diagonalization whilst results for larger systems were obtained with 
a X = 6 MERA. As L increases, the magnetizations are seen to converge toward their 
thermodynamic limit values. Results for L = 54 could not be visually distinguished from 
results for L = 18 and have been omitted in the plot. As it is characteristic of a second 
order phase transition, for large L both magnetizations develop a discontinuity in their 
derivative, with {a^ (the order parameter) suddenly dropping to zero at the quantum 



critical point (see Fig. 7.3). 



The technical details on how to numerically optimize the disentanglers and isometrics of 
the MERA to approximate the ground state |\E'o) of -f^o are analogous to those discussed 
in Chapter |5] for a \D lattice and will not be repeated here. Instead, we focus on the key 
aspect that makes the present ID scheme much more efficient than the previous 4-to-l 
scheme. For this purpose, we consider an operator Oq whose support is contained within 
a block of 2 X 2 sites of lattice £o- Direct inspection shows that, no matter where this 



block is placed with respect to the disentanglers and isometrics of Fig. |7.1[ the support 
of the resulting coarse-grained operator 0\ is also contained within a block of 2 x 2 sites 
of £i, and the same holds for any subsequent coarse-graining. This is in sharp contrast 
with the ID scheme of Refs. [EVlObl lEVlOal ICDROSj . where the minimal stable support 
of local observables (or 'width' of past causal cones) corresponded to blocks of 3 x 3 sites. 
In the present case, much smaller objects (operators acting on 4 sites instead of 9 sites) 
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are manipulated during the calculations, resulting in the announced dramatic drop in 
simulation costs. 
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Figure 7.3: Magnetizations (cr^) and (cr^) as a function of the applied magnetic field A 
for different values of the refinement parameter x- Left: Spontaneous magnetization (a^) 
for L = 54. Data fits of the form {ax) ~ (A — Xc)'^" near the critical point give a critical 
magnetic field Ac = {3.13,3.09,3.075} and critical exponent /3c = {0.320,0.321,0.323} 
for X = {2,4,6}. Current Monte Carlo estimates are Ac = 3.044 and (3^ = 0.326 |RK99l 
IBD02j . Thus accuracy increases with x- Right: Transverse magnetization (a^) for L = 6. 
TTN results for large x ^-^e taken as the exact solution (see Fig. 7.5). Whilst a x = 2 
MERA produces significantly different values, results for x = 3 ^ire already very similar 
and those for x = 6 MERA agree with the TTN solution on at least 3 significant digits. 



7.3 Results for the 2D Ising model 

We have tested the proposed scheme by investigating low energy properties of the quantum 
Ising model with transverse magnetic field, 

i^:.n.= $^4'^Wr + A5^aW, (7.4) 

(r,r') r 

on a square lattice with periodic boundary conditions (local dimension d = 2). First of 
all, we consider a sequence of lattices with increasing linear size L = {6,9, 18,54}. For 
each of them, a MERA approximation to the ground state of //ising for different values 
A G [0, 5] of the transverse magnetic field is obtained using x = 6. Computing the ground 
state for L = 54 and critical transverse magnetic field takes ~ 4 days on a 3GHz dual-core 
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desktop PC with 8Gb RAM when starting from a randomly initialized MERA|^ Fig. 7.2 



displays the expected value of the parallel and transverse magnetizations, both of which 
show characteristic signs of a second order phase transition as L increases. We emphasize 
that since the simulation costs grow only as the logarithm of L, it is straightforward to 
increase the system size until e.g. finite size effects become negligible on local observables. 



Fig. 7.3 shows how the parallel and transverse magnetizations change with increasing x, 
for L = 54. Since the cost of the simulations grows as 0{x^^), only small values of x 
can be considered in practice. However, with x = 6 one already obtains estimates for the 
location of the critical point and the critical exponent /3 that already fall within 1% of 
the best Monte Carlo results |RK99l [BD02] . 

By using the MERA to represent a two-dimensional subspace and minimizing the expec- 



tation value of -falsing, we obtain the system's energy gap AE. Fig. 7.4 shows AE as a func- 
tion of the transverse magnetic field and system size. Notice that at the critical point the 
gap closes with the system size as 1/L (dynamic exponent z = 1). Two-point correlators 



can also be extracted. Fig. 7.4 shows the correlator {ax^ax^)c = (o'lT'ai'' ') — (crx')((Jx 
along a row or column of the lattice, obtained using the scale invariant algorithm of 
Chapter |6l which directly addresses an infinite lattice at the critical point. 



7.4 Role of disentanglers 

In order to highlight the importance of disentanglers, we have also performed simulations 
with a tree tensor network (TTN) |TEV09j . This corresponds to a more orthodox real- 



space RG approach where the block of 3 x 3 sites in Fig. 7.1 is directly mapped into an 
effective site without the use of disentanglers. Recall that a 2D ground state typically 
displays a boundary law. Si ~ /, for the entanglement entropy Si of a block of / x / sites. 
To reproduce this boundary law with a TTN, one needs to increase the dimension x 
each step of the coarse-graining. Specifically, Xttn must grow doubly exponentially with 
the linear size L of the lattice. On the other hand, the cost of manipulating a 2D TTN 
grows only as a small power of Xttn- As a result, much larger values of x can be used 
with a TTN, leading to a very competitive approach for small lattice sizes |TEV09j . Fig 



7.5| (i and ii) compares the performance of the MERA and the TTN in lattices of size 



^Calculations for x = 6 are achieved by using a disentangler u with x = 4 on selected indices. The 
computation time is reduced to a few hours per point by re-using a MERA previously converged (for a 
similar magnetic field) as the starting point of a simulation. 
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Figure 7.4: Top: The energy gap as a function of the transverse magnetic field A, 
computed by exact diagonahzation for small system sizes L = {2, 3,4} and with a x = 6 
MERA for L = {6, 9}. The gap scales as 1/L at the critical magnetic field. Bottom: Two- 
point correlators {crx^crx^)c at criticality and for different values of x- The scale invariant 
MERA produces correlators that decay polynomially with the distance s = \r — r'\. As 
X increases their asymptotic scaling approaches 1/s^"'"'' with t] = 0.03 ± 0.01 |PV02] . 
Correlators have been computed at distances s = 3'^ for k = 0,1,2, .. ., where they can be 
evaluated with cost 0{x^^)- For comparison, we have included correlators obtained with 
a -D = 2 and D = 3 iPEPS |JOV"'"08] . The latter are very accurate for s = 1, 2 but decay 
exponentially after a few sites. 
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Figure 7.5: Energy error as a function of the refinement parameter x for finite systems 
of different sizes and for infinite systems. In absence of an exact solution for ground 
state energies, the errors are defined relative to the results obtained with (i) a % = 60 
TTN, (ii) a X = 9 MERA, (iii,iv) a D = 3 iPEPS |.T0V+n8j . For finite systems (i,ii), the 
MERA is compared against the TTN. The double x-axes for Xmera and Xttn have been 
adjusted so that they roughly correspond to the same computational cost. For L = 6 
the TTN is more efficient whilst for L = 9 the MERA already gives significantly better 
results. Comparison between MERA and iPEPS results for (iii) an infinite system off 
criticality and (iv) an infinite system at criticality shows very similar accuracy between 
X = 3 MERA and D = 2 iPEPS, whereas D = 3 iPEPS gives a lower (better) energy 
than X = 6 MERA. 



6x6 and 9x9. It shows that a TTN is more efficient than the MERA in computing 
the ground state of the 6x6 lattice; however, this trend is already reversed in the 9x9 
lattice, where the cumulative benefit of using disentanglers clearly outweighs the large 
cost they incur. Disentanglers, by acting on the boundary of a block, readily reproduce 
the entropic boundary law (for any value of x) and allow us to consider arbitrarily large 



systems. Fig. 7.5 (iii and iv) shows results for an infinite lattice near and at criticality. 



To summarize, we have proposed an entanglement renormalization scheme for the square 
lattice and demonstrated its scalability by addressing the quantum Ising model on systems 
of linear size L = {6,9,18,54}, with cost O(x^^logL), and on an infinite system at 
criticality, with cost 0{x^^)- The key of the present approach is the use of two types of 
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disentanglers that remove short-range entanglement residing near the corners and near 
the boundaries of the blocks while leading to narrow causal cones of 2 x 2 sites. Similar 
schemes can be built for other lattice geometries, as is demonstrated for a Kagome lattice 
in the next Chapter. 



Chapter 8 

The ground state of the spin-^ Kagome 
lattice Heisenberg model with Entanglement 
Renormalization 



8.1 Introduction 

Low dimensional spin-| quantum systems have long been the focus of intense research 
efforts, largely fueled by the search for exotic states of matter. An important example of 
a geometrically frustrated quantum antiferromagnet [Lhu05] is the spin-^ kagome-lattice 
Heisenberg model (KLHM). Despite a long history of study, the nature of its ground 
state remains an open question. Leading proposals include valence bond crystal (VBC) 
[MZOTl [SMn2l [NSnal iBAnil ISHOTI ISHnS] and spin liquid (SL) jS^^ lLE98l IWVn6l 
[Mil98l IHasOOl IHSF051 IRHLW071 IHRLW081 IJWSOSj ground states. Interest has been 
further stimulated by recent experimental work on Herbertsmithite ZnCu3(OH)6Cl2, a 
possible physical realization of the model |SNBN05] . 

Progress in our understanding of the KLHM has been hindered, as with many other 
models of frustrated antiferromagnets, by the inapplicability of quantum Monte Carlo 
methods due to the negative sign problem. Nevertheless, systems with up to 36 sites 
have been addressed with exact diagonalization |LE93[ ILLOQj , whereas the density matrix 
renormalization group (DMRG) has been used to explore lattices of order ^ 100 sites 
|JWS08] . Unfortunately, it is very difficult to infer the nature of the ground state of an 
infinite system from these results. The reason is that these lattices are still relatively 
small given the 36-site unit cell of the leading VBC proposal, or the algebraic decay of 
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correlations in some SL proposals. In larger systems, support for a SL ground state has 
also been obtained with a SL ansatz |HSF05t IRHLWOTf IHRLW08] , whereas evidence for a 
VBC has been obtained for an infinite lattice with a series expansion around an arbitrary 
dimer covering |SH07l [5H08j . but both approaches are clearly biased. 




(vii) disentangler: // disentangler: v isometryni' 



5 i| 



Figure 8.1: Coarse-graining transformation that maps (i) a kagome lattice Lq of iV sites 
into (vi) a coarser lattice L\ of iV/36 sites, (ii) The lattice is first partitioned into blocks 
of 36 sites, (iii) Disentanglers u are applied across the corners of three blocks, followed 
by (iv) disentanglers v applied across the sides of two neighboring blocks, (v) Isometrics 
w map blocks to an effective site of the coarse-grained lattice, (vii) Tensors m, v and w 



have a varying number of incoming and outgoing indices according to Eq. 8.1 



In this Chapter we report new, independent numerical evidence in favor of a VBC ground 
state for the KLHM model obtained by a numerical study with entanglement renormal- 
ization. Entanglement renormalization is a real space RG approach which, through the 
proper removal of short-range entanglement, is capable of providing an approximation to 
ground states of large ID lattices, as has been demonstrated in Chapters |3| |4] and [7], by 
means of a multi-scale entanglement renormalization ansatz (MERA). After describing 
a MERA scheme for the Kagome lattice with periodic boundary conditions, we address 
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lattices of = {36, 144, oo} sites. Our simulations converge to a VBC state compati- 
ble with that first proposed by Marston and Zeng |MZ91j and revisited by Nikolic and 
Senthil |NS03] . and by Singh and Huse |SH07l ISHOSj . For an infinite lattice we obtain 
an energy per site E=-0.4322. This energy corresponds to an explicit (MERA) wave- 
function and therefore provides us with a strict upper bound for the true ground state 
energy. Importantly, its value is lower than the energy of any existing SL ansatz on a 
sufficiently large lattice, which we interpret as strong evidence for a VBC ground state in 
the thermodynamic limit. These results also demonstrate of the utility of entanglement 
renormalization to study 2D lattice models that are beyond the reach of quantum Monte 
Carlo techniques. 



8.2 ER on the kagome lattice 



The present approach is based on the coarse-graining transformation of Fig. |8.1[ which is 
applied to a kagome lattice Cq made of N sites. It maps blocks of 36 sites of Cq onto single 
sites of a coarser lattice Ci made of A^/36 sites. A Hamiltonian Hq defined on lattice Cq 
becomes an effective Hamiltonian Hi on lattice Ci. Analogously, the ground state |\I^o) 
of Hq is transformed into the ground state oi Hi. The transformation decomposes 
into three steps. Firstly disentanglers u, unitary tensors that act on 9 sites, are applied 
across the corners of three neighboring blocks. Then disentanglers v are applied across 
the sides of two neighboring blocks; these tensors reduce ten sites (each described by a 
vector space C2 of dimension 2) into two effective sites (each described by a vector space 
of dimension x). Finally isometries w map the remaining sites of each block into a 
single effective site of Ci. Thus the tensors u, v and w, 
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v'^v = 












w'^w = 





(8.1) 

transform the ground state |\l/o) of lattice £0 into the ground state of lattice £1 
through the sequence 

|vl/o) A|M/'o) A|M/'o')^|M/i). (8.2) 

The disentanglers u and v aim at removing short-range entanglement across the bound- 
aries of the blocks; therefore states |\E'q) and |\1/q) possess decreasing amounts of short 
range entanglement. If state |\l/o) only has short-range entanglement to begin with, then 
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it is conceivable that the state has no entanglement left at all. For a finite lattice 
(A^ = 144) we consider a state |\E'o) that after the coarse graining transformation give 
rise to an entangled state on A^/36 = 4 sites. For an infinite lattice we will instead 
make an important assumption, namely that |\E'i) is a product (non-entangled) state. 
How short-ranged must the entanglement in |\E'o) be for this assumption to be valid? By 
reversing the transformation on a product state I^E'i), it can be seen that each site in |\E'o) 
is still entangled with at least 84 neighboring sites. 



Figure 8.2: The 36-site unit cell for the honeycomb VBC, strong bonds are drawn with 
thick lines. Three different types of strong bonds can be identified; the six bonds belonging 
to the pinwheels (red), six bonds belonging to each 'perfect hexagon' (blue) and the 
parallel bonds between perfect hexagons (green). Dotted arrows indicate the axis where 
spin-spin correlators have been computed. Bond-bond correlators have been computed 
between the reference bond (1) and the other numbered bonds. 



8.3 Results and discussion 

The disentanglers and isometries {u, v, w) were initialized randomly and then optimized 
so as to minimize the expected value of the KLHM Hamiltonian, 





(8.3) 



by following the algorithms described Chapter [sj with cost 0{2^'^x^x'^) Q Specifically, 
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for lattices with = 36 and 144 sites, the resulting (one-site and four-site) Hamiltonian 
Hi is diagonalized exactly. Instead, for = oo, we use the finite correlation range 



algorithm (Sect. 5.4.4). All computations led to highly dimerized wave-functions of the 



VBC type. In order to explain the results, consider the exact 'honeycomb' VBC state. 



denoted |h-VBC), whose 36-site unit cell is shown in Fig. 8.2, Each unit cell contains two 
'perfect hexagons' (resonating bonds around a hexagon) and a 'pinwheel'. Three different 
types of strong bonds can be identified: those of the pinwheels (red), parallel bonds (green) 
and perfect hexagons (blue). The pinwheel and parallel bonds are singlets (energy per 
bond = —0.75) while the perfect hexagons are in the ground state of a periodic Heisenberg 
chain of 6 sites (energy per bond = —0.4671). The rest of links have zero energy. We call 
a 'honeycomb' VBC a state that has strong bonds according to the above pattern, even 
though the rest of bonds (weak bonds) need not have zero energy. The 'honeycomb' VBC 
was originally proposed by Marston and Zeng jMZQlj (see also |NS03l [SH071 ISH0"8] ) . Our 
simulations with A^ = 144 and N = oo produce a VBC of this type as the best MERA 
approximation to the ground state. 



The energies obtained for an infinite lattice are shown in Table 8.1 For each value of 
X, the MERA is an explicit wave-function and therefore provides an upper bound to the 
exact ground state energy. Energies computed for the A^ = 144 lattice matched those 
of the infinite lattice to within 0.02% and have been omitted. These N = oo energies 
also match closely those obtained by series expansion in Ref. |SH07j . and are lower than 
those obtained in Ref. |.TWS08j with DMRG {E = -0.43160 for A^ = 108) and in Ref. 



^The computational cost of simulations can be signifcantly reduced by incorporating symmetries of 
the system into the MERA tensor network [SPV09 . In this work, a U(l) symmetry is incorporated into 
the MERA to allow large x simulations that would otherwise be unaffordable. 



Table 8.1: Ground state energies as a function of x- 



X 


A^ = oo 


A^ = 36 


A^ = 36 






(rand init) 


(|h-VBC) init) 


2 


-0.42145 


-0.42164 


-0.42143 


4 


-0.42952 


-0.42816 


-0.42715 


8 


-0.43081 


-0.43199 


-0.43148 


12 


-0.43114 
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[RHLWOT] with fermonic mean-field theory and Gutzwiller projection {E = —0.42863 for 
N = 432). We further notice that, where finite size effects are still relevant, such as in 
the = 108 case, they tend to decrease the ground state energy. 
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Figure 8.3: (top) Bond energies for the 36-site unit cell of infinite MERA wave-functions, 
for two different values of x, as compared to those of an exact honeycomb VBC, |h-VBC), 
and those of a spin liquid, which by definition has all equal strength bonds. The MERA 
wave-functions clearly match the proposed honeycomb VBC; we identify (i) the six strong 
'pinwheel' bonds (red bonds), the six 'parallel' bonds (green bonds) and (iii) the 12 'perfect 
hexagon' bonds (blue bonds). The (iv) remaining 48 bonds are the weak bonds of the unit 
cell, (bottom) Bond energies for the 36-site lattice. Here a randomly initialized MERA 
converges to a dimerized state that does not match the honeycomb VBC pattern, but 
gives lower overall energy than a honeycomb VBC initialized MERA of the same x- 



Fig. 8^ shows the distribution of bond energies obtained for the N = oo lattice. With 
X = 4, one observes an energy increase per site over |h-VBC) of ~ 0.08 in the parallel 
(green) bonds and also in some of the hexagon (blue) bonds, with the weak bonds having 
lower energy in return. As x is increased, the energy of the 'strong' bonds becomes slightly 
larger and that of the 'weak' bonds continues to decrease. However, the dimerization 
clearly survives: the bond energies are not seen to converge to a uniform distribution as 
required for a SL. 
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Figure 8.4: Spin-Spin correlators along arrows 'A' and 'B' of Fig. |8.2| for infinite lattice 
MERA of X = 4 and x = 16. Although along both lattice directions considered the 
correlators decay exponentially, the decay along arrow 'A' (a line joining two perfect 
hexagons) is seen to be slower than along arrow 'B' (a line joining perfect hexagon to 
pinwheel). The plateaus marked (i), (ii) and (iii) show the correlation is the same with 
both spins of a strong bond. 



Fig. 8.4 shows spin-spin correlators evaluated along two different lattice axis A and B (cf. 
Fig. 8.2) for = oo. These correlators decay exponentially with well defined 'plateaus', 
where the correlation is the same with both spins of a strong bond. Correlations along 
the line joining a perfect hexagon and a pinwheel are seen to decay faster than along the 



line joining two perfect hexagons, consistent with the observation from Fig. 8.3 that the 
perfect hexagon bonds remain almost exact singlets even for high values of x- Table 8.2 
shows bond-bond connected and disconnected correlators, Cia and Dia, 



. {[S-S)^[S-S)J (8.4) 
= Ci,^~((s-S))((s-S) ), (8.5) 



between a reference bond '1' and a surrounding bond a = 1, ■ ■ ■ ,14 (cf. Fig, 8.2). While 
disconnected correlators decay exponentially with distance, some connected correlators 
remain significant at arbitrary distances, demonstrating the long-range order of the VBC 
state. 



Let us discuss the results for a lattice with A^ = 36 sites. When initialized with random 
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Table 8.2: Bond-Bond correlators for x = 16 {N = oo). 
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tensors, the MERA produced VBC type configurations which typically did not match 
the honeycomb VBC, although simulations initialized in the state |h-VBC) retained the 
honeycomb VBC pattern (see Fig. 8.3). Here the randomly initialized VBC produced 
a lower energy (0.5% above the exact diagonalization result E = —0.438377 of |LE93j ) 
than the honeycomb VBC type solution for an equivalent value of x (cf- Table 8.1). 
These results strongly suggest that finite size effects in the = 36 site lattice lead to a 
significant departure from the physics of the infinite system. 



8.4 Defective valence bond crystals 

In this section we detail a modification of the MERA optimization algorithm of Sect. 



5.4| that was found necessary to ensure proper convergence of some simulations of the 
KLHM. All simulations of the KLHM we found to converge to highly dimerized VBC 
type states. For lattices of A^ = {144, oo} sites the lowest energy state that was found, 
or best approximation to the ground-state, had a pattern of strong bonds matching the 
'honeycomb' VBC state. However, not all randomly initialized simulations converged to 
this state; in some instances the MERA wave-function converged to a 'defective' VBC 



state, an example of which is shown Fig. 8.5 The observed defective VBC's ranged from 



states which matched the honeycomb VBC with only a few misplaced strong bonds to 
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states with an almost completely disordered placement of strong bonds. A defective VBC 
was observed to have on average 0.05% higher energy per site than the corresponding 
honeycomb VBC type MERA of the same x- The exact difference in energy depended on 
the amount of defects; typically, states that deviated more from the honeycomb VBC had 
higher energy than states which deviated less. In the instances that the method did not 
converge to the honeycomb VBC one may conclude that the optimization of the MERA 
became trapped in a local minimum. 

Becoming trapped in a local minimum is not surprising given that the optimization of the 
MERA is based upon updates of individual tensors. During the initial iterations of the 
optimization, locally stable structures, such as singlets between neighboring spins, form. 
But they do not necessarily form in places compatible with a global minimization of the 
energy. Once the method has converged to a particular VBC state, the transformation 
required to bring the state to a different VBC, with a new pattern of strong bonds, requires 
simultaneous shifts of many strong bonds. This is unlikely to occur with an optimization 
based on energy-lowering updates of individual tensors. 

There are many ways to decrease the risk of becoming trapped in a local minimum due 
to the formation of local singlets. In the present work this is achieved by restricting 
how much a tensor is allowed to change in one single update; this may be thought of as 
decreasing the rate at which the MERA is 'cooled-down' from an initial (high-energy) 
state to the (low-energy) approximation to the ground-state. The goal of this restriction 
is to prevent locally stable structures from forming too early in the optimization. 

In order to discuss how this restriction may be implemented, we first refresh some details 
of the variational MERA optimization, see Sect. 5.4[ In the standard algorithm, a tensor 



w in the MERA is replaced with a better (that is leading to lower energy) tensor w' as 
follows. First one computes the linearized environment of w, denoted T^. Then, given the 
singular value decomposition of the environment, T^^, = USV'', one chooses the updated 
tensor to be w' = VW . Here we wish to restrict the amount e the tensor can change, 
\\w — w'W < e. This is achieved by using a modified environment defined as 

t„ = + \w (8.6) 

for some A > 0. As before, given the singular value decomposition of the modified 
environment T^^ = U SV'^ , the updated tensor w' is chosen to be w' = VW . The parameter 
A acts as a soft constraint to ensure the tensor change e is kept small with each update, 
a larger value of A giving a smaller change e. In the limit of A — )■ oo, the tensor w would 
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remain constant, e = 0, in the update. Simulations of the KLHM with a suitable value of 
A, say A > 10 for the initial optimization, were much more likely to converge to a defect 
free honeycomb VBC than were unrestrained (A = 0) optimizations. 



N = 144, defective valence bond cystal 



N = 144, honeycomb valence bond cystal 





Figure 8.5: Converged x = 8 MERA solutions for the N = 144 site KLHM. The thickness 
of a line between two lattice sites is proportional to the absolute bond energy, \E\. Blue 
solid lines are used for negative energy bonds, E < 0, and red dashed lines for positive 
energy bonds, E > 0. (Top) An example of a 'defective' VBC. Local structures that 
match the 'pinwheels' and 'perfect hexagons' of the honeycomb VBC are observed, but 
long range order is absent. (Bottom) By restricting the rate of lattice 'cooling', as per 



Eq,8^, one more reliably obtains the honeycomb VBC from a randomly initialized MERA. 
On average a defective VBC had 0.05% greater energy than the honeycomb VBC of the 
same value of x- 



8.5 Conclusions 

To summarize, we have used entanglement renormalization techniques to obtain new 
numerical evidence indicating that the ground state of the KLHM is of the honeycomb 
VBC type. In order to assess the robustness of this result, we briefly discuss some of the 
limitations of the present approach. 



Firstly, the coarse-graining transformation of Fig. 8.1, which maps 36 sites into one site, 
was designed to ensure compatibility with the 36-site unit cell of honeycomb VBC type 
solutions. While our approach did not preclude solutions with a smaller, compatible unit 
cell (such as a 12-site unit cell or a fully translation invariant solution), we cannot rule 
out the possibility that a state with an incompatible unit cell might have a lower energy. 
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Secondly, the infinite lattice was investigated by restricting the range of entanglement in 



the ansatz to blocks of 84 spins, imposed through an unentangled state in Eq. 8.2 
This restriction was only implemented after preliminary simulations with x = 12 had 
produced identical energies irrespectively of whether was allowed to be entangled. 
However, it could still be that entanglement in |\l'i) would make a big difference for larger 
values of x- We find this scenario quite unlikely, but could not test it due to computational 
limitations. 



Finally, the MERA is an essentially unbiased method provided that the candidates to be 
the ground state of the system have all a relatively small amount of entanglement. But 
when deciding between a VBC (which mostly has short-range entanglement) and e.g. the 
algebraic SL of Refs. jRHLWOTf IHRLW08"] (significantly more entangled at all length 
scales), it might well be that the MERA is biased toward the low entanglement solution. 
Therefore our results do not conclusively exclude a SL ground state. We emphasized, 
however, that the ground state energies obtained with the MERA are lower than the SL 
energies of Refs. [RHLWOTt lHRLWn8l iJWSnS] . 



Future work includes the computation of singlet and triplet excitation gaps, and studying 
the effect of adding a magnetic field and anisotropic terms to the Hamiltonian. Such 
additions to the Hamiltonian can be considered without modifying the present approach. 
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Chapter 9 

Conclusion and Future Outlook 



9.1 Summary 

The content and major results of this thesis are summarised below. Firstly, Chapter [2]re- 
viewed the foundations of entanglement renormalization and the MERA as first introduced 
in Refs. |Vid07tlVid08j . Chapters [s] and |4] then explored the performance of entanglement 
renormalization in the simple setting of free fermions and free bosons respectivly. Here 
it was demonstrated that a real-space renormalization group transformation based upon 
entanglement renormalization could induce a sustainable coarse-graining transformation 
in variety of D = 1, 2 dimensional free fermion and free boson models. The substainability 
of the RG transformation with entanglement renormalization allowed investigation of ar- 
bitrarily large D = 1,2 dimensional lattice systems without loss of accuracy, as was shown 
through comparison with exact solutions. An exception was found for critical 2D free- 
fermion models with a ID fermi surface, whose entropy is known to scale as Sl ~ Llog(L), 
where ER could no longer produce a be sustainable RG transformation. It was concluded 
that a generalized MERA would be required to analyse such systems. 

Chapter [5] described how to compute expected values of observables from a MERA and 
also presented optimization algorithms to compute the ground-state MERA for an arbi- 
trary local Hamiltonian. A highlight of algorithms was their computational cost which, 
for a chain of sites, scales as 0{N) for non-translation invariant systems, 0(log(A^)) 
for translation invariant systems and scales independant of for scale-invariant systems. 
This Chapter also presented benchmark results for several ID spin models, mainly for 
Ising and Potts models, including accurate computation of ground energy, local observ- 
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ables, correlators and energy gaps. Chapter |6] analysed the connection between entangle- 
ment renormalization and conformal field theory (CFT). It was demonstrated explicitly 
for several ID spin models that this connection could be exploited to extract most of 
the conformal data of the CFT that described the models in the continuum limit. This 
included accurate calculation of the scaling dimensions, the central charge and OPE co- 
efficients for Ising and Potts models. 

In Chapter [7| the optimisation algorithms presented in Chapter |5] were generalised to the 
specific problem of analysing quantum many-body problems on D = 2 dimensional lat- 
tices. Benchmark calculations with the quantum Ising model demonstrated the ability of 
the MERA to give accurate ground state properties on arbitrarily large (or infinite) lattice 
systems for a non-integrable model in D = 2 spatial dimensions. The benchmark results 
included accurate calculations of critical exponents, expected values of local observables, 
two-points correlators and energy gaps. Finally, Chapter [s] applied ER to study the spin-| 
kagome lattice Heisenberg model (KLHM), a long-standing problem in condensed matter 
physics. The results of this Chapter, which is the first demonstration of ER to study an 
open problem in condensed matter physics and a problem beyond the reach of quantum 
Monte-Carlo, offers convincing evidence that the ground-state of the infinite KLHM is a 
honeycomb valence bond crystal. In particular, the energy of the MERA approximation 
to the ground-state of the infinite system, which serves as an exact upper bound to the 
true ground energy, is significantly better than any such bound previously produced. 



9.2 Future Work and Outlook 

Entanglement renormalization offers a non-perturbative means to investigate quantum 
many-body systems, as such ER has a diverse range of applications that remain to been 
explored or further developed. 

Firstly one could consider extensions to the research of Chapter |6] (see also |GMF08l 
IPEVOQt IMRGF09] ) which made a connection between the MERA and conformal field 
theory (CFT). The procedure to compute the local scaling operators (and their corre- 
sponding scaling dimensions) of a critical theory, as described Chapter [6| could be gen- 
eralised to obtain a class of non-local scaling operators. Thus it might be possible to 
obtain a complete set of the scaling operators which characterise a CFT. It also remains 
to generalise the ER algorithm to the study of types of conformally invariant systems 
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other than the ID infinite chains previously considered. These may include systems with 
a boundary, homogeneous systems with a defect (or several defects), an interface between 
two different critical systems, or a Y-intersection of three spin chains. 

There is much potential for the use of ER to study 2D quantum systems. In particular 
ER may have a big impact in the areas of frustated magnets, systems of interacting 
fermions and topologically ordered systems, as these are problems for which previously 
there has often not been any satisfactory investigative techniques. The results in Chapter 
[8] for the kagome lattice Heisenberg model are a demonstration of the potential of the 
method for frustrated magnets. Entanglement renormalization could be similarly applied 
to investigate many other models of magnetism that remain not well understood; these 
may include Jl — J2 models, ring exchange models and quantum dimer models. In 
particular, it is possible that ER could be used to study the existance of a spin-liquid phase 
in a 2D magnetic system, a problem that is of much interest to the condensed matter 
community. Recent developments have led to a generalisation of the ER optimisation 
algorithm to allow the investigation of systems of interacting fermions |CEW10t ICV091 
IPBElOt IBPEOQj . At present most of the studies of fermions with ER have been tests to 
benchmark the algorithm, and there remains a wide range of interesting fermionic models 
that could potentially be investigated with ER. In regards to topological order, previous 
studies also have shown that ER can be used to characterise the infra-red fixed point 
of topologically ordered systems |AV08i IKRVOQj . Further research could utilise ER to 
investigate, for instance, the stability of a topological phase under perturbation. 

Prior to the further application of ER to study 2D lattice systems it may be necessary 
to devise ways to improve the computational efficiency of the algorithm, which typically 
scales as a large power of the bond dimension x- Significant reductions in computational 
cost can come from incorporation of global or local symmetries into the tensor network; 
research into the incorporation of global symmetries is already well progressed |SPV09] . It 
may also be possible to reduce costs by introducing approximations into contraction of the 
MERA tensor-network or to incorporate monte-carlo type sampling into the algorithm. 

In the long term, it is possible that entanglement renormalization, as a tool for analysing 
many-body systems, may find applications beyond the more immediate applications in 
condensed matter physics. Applications could arise wherever quantum many-body effects 
are important; for instance in areas such as string-theory, quantum gravity and quantum 
chromo dynamics. Entanglement renormalization allows the evaluation of the continuum 
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limit of a lattice model, hence may be useful for analysing problems that would be difficult 
to analyze in the continuum directly. Recent work on the holographic principle |Swi09j 
might be a first step in this direction. 



Appendix A 

The ground-state of free fermion models in 
two dimensions 



A.l Particle conserving free-fermions 

In this appendix we derive an expression for the ground-state correlation matrix of free- 
fermion models on a 2D square lattice (a thorough derivation of ground-state correlation 
matrices for ID free-fermion models can be found e.g. |LRV04] ) . The ground-state 
covariance matrix is the starting point for the analysis of free fermions with entanglement 
renormalization considered in Chapter |3| We begin with by considering nearest neighbour 
free-fermion models that have only particle conserving terms on a 2D lattice of N x N 
sites 



ri,r2=-(Af-l)/2 

We now perform a 2D Fourier transform of the fermionic operators a into a new set of 
operators h defined by 
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The inverse transformation is similarly defined 

1 
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(A.3) 



By exploiting the orthogonality properties of the (unitary) fourier transform, it is seen 



that the hopping terms in Eq. A.l become on-site interaction terms in the fourier basis 
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(A.4) 
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The on-site interaction term in Eq. A.l which may be interpreted as a chemical potential, 
remains invariant under the fourier transform 

^ri,r2^ri,r2 = ^^kiM^kiM' (^-5) 

ri,r2 ki,k2 

On subsitution of Eqns. A.4 and A. 5 the original Hamiltonian of Eq. |A.l is recast into 
diagonal form 

{N-l)/2 

H= ^k^M^l^M^k^M (A.6) 

fci,fc2=-(Af-l)/2 

with disperion relation A given as 

Afc^,fc2 = - cos(27rA;i/iV) - cos(27rA;2/iV) - A. (A.7) 

The ground-state of Hamiltonian |A.6 is defined as the state with all negative energy 
modes, ^k^M < 0? occupied (the fermi sea) and the rest of the fermionic modes empty. It 
follows that the ground-state correlation matrix, in the fourier basis, is given by 
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(A.8) 



The inverse transformation relations of Eq, A.3 can be used to transform the correlation 
matrix into the original (spatial) basis 
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where T denotes the fermi sea. In order to take the thermodynamic hmit of Eq,A.9, we 
first define new variables 

where, by definition, 0i,02 ^ [~7r,7r]. The thermodynamic limit, — )■ oo, is now taken 

<o«n,r.)„ = ^ [ e''-'^^e''''^'dM<p2 (A.ll) 



47r" 



where the integral is taken over those modes within the fermi sea J-". Explicitly, the integral 
should be evaluated over the region A > 0, where A^^ = — cos(0i) — cos(02) — A. In 
most cases the intergal in Eq jA.ll| cannot be evaluated analytically and one must resort to 
numerical integration to obtain approximate values for the correlators. However, for the 
case of half-filling (A = 0), the fermi-sea has a particularly simple form and the integral 
may be carried out analytically as follows. Firstly, a change of variables is made 

= 4>1 + 4>2, V2 = (t>l-4>2- (A. 12) 



Expressed in the new variables, the integral of Eg A.ll can be evaluated 
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A. 2 Non particle conserving free-fermions 



We now consider more general quadratic models of fermions than was initially considered 
with Hamiltonian A. 1 Specifically, we consider Hamiltonians with non particle conserving 
terms of the form aa and a^a^ . The generic nearest neighbour Hamiltonian on an A^ x A^ 
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lattice is written 
H = 
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(A.15) 

This model of free-fermions depends on a parameter A, with represents the chemical 
potential, and another parameter 7, which represents the anisotrophy. The free-fermion 



model of Eq. A.l corresponds to setting the anisotrophy term to zero in the present 
model. Similar to the derivation of last section, we shall first proceed with a fourier 
transform of the fermionic operators a. However, as we shall see, the Hamiltonian of 



Eq. A.15 also requires an additional Bogoliubov transformation to bring it into diagonal 



form. Applying the Fourier transform of the fermionic operators defined Eq. |A.2[ the non 



particle conserving terms of Eq. A.lS transform as 
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while the other terms in the Hamiltonian transform in the same manner as Eqns. A. 4 



and A. 5 The Hamiltonian expressed in the Fourier basis b is given 
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While the Fourier transform has eliminated most of the off-diagonal couplings from the 
Hamiltonian, couplings of the form bk^^k2b-ki,-k2 have been introduced. We now require 



a Bogoliubov type transformation to fully diagonalise the Hamiltonian of Eq. A.17 The 
Bogoliubov transformation, which acts independantly for every value {ki, /C2), involves tak- 
ing linear combinations of the annihilation/creation operators 6^, b to form new fermionic 
operators c^, c and is defined 
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with 6'fci,fc2 as a yet undefined function of {ki, k2). In order to determine tlie function Oki^k2 
it is first assumed tliat tliat Hamiltonian is diagonal in the Bogoliubov basis c 

(N-l)/2 
kr,k2=-(N-l)/2 

and then one works backwards to match the coefficients of Eq. A.20 with those of Eq. 



A. 17 Using Eq. A. 19 
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On subsitution of this expression into the Hamiltonian of Eq. A.20 we get 
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of the trigonometric identities 



(A.22) 



where we have made use 

cos(x) = cos^(x/2) — sin^(x/2) 
sin(x) = 2 sin(x/2) cos(x/2). 

By matching the coefficients of Eq. A. 17 with Eq. |A.22 
the transform weights is 



is seen 



cos (6',fcj^fc2) = ^ki,k2/^ki,k2 
sin (6'jtj^fc2) = ^ki,k2/^ki,k2 

where the dispersion relation fl may be written as 



(A.23) 

that the correct choice for 
(A.24) 



'fcl,A:2 



The Hamiltonian of Eq. A. 15 has now been tranformed into a diagonal form 



{N-l)/2 



H- ^ ^kuk2Ck^,k2'^ki,k2- 

ki,k2=-(N-l)/2 



(A.25) 



(A.26) 
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Note that, due to the action of the Bogohubov transformation which mixes creation and 



annihilation operators, the dispersion in Eq. |A.26| is positive defined, unhke that of Eq 



A. 6 from the previous section. Now that the correct transforms to diagonahse the free- 
fermion Hamiltonian has been identified, we turn our attention to finding the ground-state 
correlators. The ground-state correlation matrix is defined in the c basis as 

From this definition we can derive expressions for the ground-state correlators ('aJo'^' 
in the original (spatial) basis 

{N-l)/2 

'a) a \ - — V (h^ bu u \e'^^iiKri+k'2r2) 

ki,k2,k['k'2=-{N-l)/2 



00"'nr-2 , 



^ Yl ((^sin(^fc,,fc2/2)c-fci,-fc2 +cos(0fc,,fc2/2)4^ 

ki,k2,k'^,k'2 

cos(4.,fc,/2)cfc.,fc, - ^sin(^,,,fc,/2)cl,, e2-('=i'^i+'=2'-2) 



k\,k2 



k\,k2 



The other correlation matrix (aLaJ ) is similarly derived 

\ ^ V GS 

flt t \ _ J_ V Ih'^ h,, ,,\p-^^i(k[n+k',r2) 
a00«rir2/ - \'kiM'^k[,k'^ J ^ 



5Z ((^sin(0fc,,fej2)c_fc,,_fc2 +cos(^fc,,fe2/2)4^fcJ 

(isin(^,,,fc,/2)c_fc.,_,. + cos(^^,,./2)4, e-2-^«^^+'=^^^) 

^ ^sin(4„,,/2)cos(^.,.,/2) (c-fe,-.,4, ,,)e-2-(^'i^^+^^'-^) 

^ ^ sin(4„,,/2) cos(^,„,,/2)e2'^^('=^^^+^^'^^) 

fci,/C2 

^ ^ / ^^^'^(^fci.fc2) \ 27r»(fciri+fc2r-2)_ (A. 29) 

^ fci,fc2 V ^ / 
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Finally, for clarity, we restate the Fourier/Bogoliubov transform weights 



COS 



sin(6'fc, 



k2, 



^ki,k2 



(A.30) 



with A defined Eq. A. 7 and $ defined Eq. A.18[ As with the previous section one could 
take the thermodynamic limit to obtain integral equations for A.28 and A.29 However, 
as the resultant integrals do not have an analytic solution, it is convenient to leave these 
equations as summations. 
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Appendix B 

The ground-state of harmonic lattice 
systems 



In this appendix we diagonalise a nearest-neighbour harmonic lattice in order to derive 
an expression for the ground-state covariance matrices of the system. Similar derivations 
may also be found e.g. |AEPW02| IPEDC05| ISkrOSj . The ground-state covariance matrix 
is the starting point for the real-space RG analysis considered in Sect. 4.4.2 of Chapter 
|4j We also offer a procedure for regularising the case of zero-mass, which would otherwise 
be divergent. As introduced in Eq. |4.1[ we focus on Harmonic systems with Hamiltonian 

(B.l) 



H 



(N-l)/2 
=-{N-l)/2 



Recall that operators Pr and are the usual canonical coordinates with commutation 
[Prji^r'] = ihSrr'- By a Fouricr transform of the canonical coordinates 



r=l 
1 ^ 

... 



-2TTirK/N 



-2TTirK,/N 



r=l 



the Hamiltonian of Eq. B.l is bought into diagonal form 



H 



{N-l)/2 

E 

=-{N-l)/2 



Pi 



mV + 8irsin^(7rK/A^) g, 



(B.2) 



(B.3) 



It is seen that Eq. B.3 represents a system of uncoupled oscillators, each with an 
effective mass dependant on k. It follows that the ground-state in this basis is the product 
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of the single oscillator ground-states 

1 



(PreiPK2)GS 



1 



(B.4) 



'm2w2 + 8K sin2(7rKi/A^) 
The covariance matrices in the original (spatial) basis are derived by using the inverse of 



the Fourier transforms of Eq. |B.2| to obtain 

' m^Lo^ + 8K sin^ { —) 



2 



{PoPr) 



imr) 



GS 



OS 



1 

2iV S 



cos 



1-JV 
'~ 2 

2 



2Ar ^ 



cos (?f^) 



vat; 



l-JV 
2 



(B.5) 



^m2w2 + sin^ (^) 

As we are interested in the systems of infinite extent, it is possible to take the thermo- 



dynamic (A^ — )■ oo) limit of Eq. B.5, in which the sums will be replaced by an integrals. 
However, in all but a particular case, to be addressed shortly, the resulting integral equa- 
tions cannot be solved analytically. It is often more convenient, if one desires correlators 
from the infinite system between sites at most R sites distant, to use the finite A^ equa- 



tions |B.5| with A^ ^ i? to in order to compute approximate correlators for the infinite 
system. For the m = case it is possible to evaluate the integral equations exactly; taking 
the thermodynamic limit of Eq. 



B.5 



(PoPr) 



with k 

— TT 

2K 



Ottk 

~N~ gives 



GS 



2ti 



COS {kr 

TV 

COS {kr 



sm 



dk 



GS 



sm 



(I) 



-dk. 



(B.6) 
(B.7) 



The integral for the p-quadrature evaluates as 



{PoPr) 



2K 



GS 



TT 



1 



2r + 1 2r - 1 



(B.8) 



Unfortunately, the integral for the g-quadrature in Eq. B.7 is divergent, {^q)Qr = oo , 
for all r. One can proceed by regularizing the integrals with a small momentum cut-off 
that is only modes with > e are integrated. The limit as £ — )■ is then investigated. 
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To evaluate the g-quadrature integral of Eq. B.7 first the substitution x = e*'^/^ is made 

1 



7rV2K 



X 



2r 



-dx — 



— 1 J x"^ — 1 

k=e k=—TT 



X 



2r 



dx 



(B.9) 



The integrand in Eq. B.9 is known to have a finite series expansion 

s=l 



X^ 



x2 - r 



(B.IO) 



Using Eq. B.IO one may spht the integral of Eq. B.9 into a convergent quantity /(r 



(which contains the correlations) and a divergent constant fi^ (which is independent of 



r). Thus the integral of Eq. B.9 is evaluated to give 
with spatial correlators /(r) defined 



fir) 



27tV2K s=i 

and the constant (which divergent in e) defined 



= ; log f cot f ^ 



(B.ll) 



(B.12) 



(B.13) 



By using the regularized expression for correlators of Eq. B.ll quantities such as the 
difference between correlators are seen to remain finite in the limit of e taken to zero and 
may be evaluated exactly 



lim ((gogr)GS;e ' (mv) GS;e) = f (r) ' f (r) . 



(B.14) 



Similarly it can be shown that, although the entanglement entropy of a block of L modes 
diverges in the massless case, the difference in entropy between two blocks is also conver- 
gent in the limit as e is taken to zero. Thus, for instance, in the zero-mass system we may 
still make sense of how the entropy of a block of L modes changes along the RG fiow as 



was demonstrated Fig. 4.8 
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